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PREFACE 


This document contains the proceedings of the Workshop on Computational Modeling of 
Tires held at the Peninsula Graduate Center of the Virginia Consortium of Engineering and 
Science Universities, October 26-27, 1994. The workshop was jointly sponsored by the 
University of Virginia Center for Computational Structures Technology and NASA. 

Workshop attendees came from government agencies, tire industry, commercial software 
companies, and universities. The objectives of the workshop were to assess the state-of- 
technology in computational modeling of tires and to provide guidelines for focused future 
research leading to an enhanced capability for the numerical simulation of tire response. 

Certain materials and products are identified in this publication in order to specify 
adequately the materials and products that were investigated in the research effort. In no case 
does such identification imply recommendation or endorsement of products by NASA, nor 
does it imply that the materials and products are the only ones or the best ones available for this 
purpose. In many cases equivalent materials and products are available and would probably 
produce equivalent results. 


Ahmed K. Noor 

Center for Computational Structures Technology 
University of Virginia 
Hampton, VA 

John A. Tanner 

NASA Langley Research Center 
Hampton, VA 
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INTRODUCTION 


A significant capability has been developed over the past two decades for the 
computational modeling and analysis of tires. Several in-house and commercial software 
systems are currently available for accurate determination of static load deflection 
characteristics, contact area shape, free vibration modes, and quasi-static rolling resistance. 

This success can be attributed, in part, to the joint NASA/U.S. tire industry National Tire 
Modeling Program (NTMP). This program was initiated in 1983 and was completed in 1991. 
At its peak seven U.S. tire companies, four universities, and a software developer participated 
in the program. The NTMP established a set of benchmark tire modeling problems, developed 
a national database of experimental tire measurements, and led to the development of a number 
of tire modeling strategies and solution algorithms that could be incorporated into new or 
existing software packages. 

Despite the significant advances made in the computational modeling and analysis of 
tires, a number of technology needs and related tasks must be addressed by the research 
community to enhance the state-of-the-art in computational modeling of tires and to make it an 
effective tool in tire design. The joint University of Virginia/NASA workshop held at the 
Peninsula Graduate Center of the Virginia Consortium of Engineering and Science 
Universities, October 26-27, 1994 focused on the status of computational modeling of tires and 
the current and future needs for further development of this technology. The list of pacing 
items given in this introduction was compiled from a number of participants. 


COMPUTATIONAL NEEDS 

The technology needs identified by the participants can be grouped into the following five 
major headings: 1) material characterization of filled rubber and heat generation mechanisms; 

2) understanding the physical phenomena associated with frictional contact, wear and noise 
generation; 3) accurate determination of tire footprint shape, size and pressure distribution; 4) 
prediction of wear resistance, failure and fatigue; and 5) tire/vehicle interaction and its effect on 
riding and han dlin g qualities and performance. 

For each of the aforementioned items, attempts should be made to exploit the major 
characteristics of high-performance computing technologies, as well as the future computing 
environment. The five primary technology needs and related tasks are described subsequently. 

1. Material characterization of cord and filled rubber and heat generation mechanisms. 
The reliability of the predictions of the tire response, damage and failure is critically dependent 
on the accurate characterization and modeling of the material behavior and heat generation 
mechanisms under different operational conditions. 

2. Understanding the physical phenomena associated with frictional contact, wear and 
noise generation. This includes understanding the frictional contact between the rubber and 
road and its effect on handling and braking forces, and mechanisms of wear and abradability of 
material. 

3. Accurate determination of the footprint shape, size and pressure distribution which 
includes development of simple and fast numerical algorithms for handling continuous contact 
and modeling the actual tread patterns. 

4. Prediction of wear resistance, failure and fatigue. This includes capability for the 
accurate calculation of contact forces, and predictions of failure/fatigue and wear initiation and 
propagation based on stress/strain history. 
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5. Tire/vehicle interaction and its effect on riding quality and handling performance - 
relating the response characteristics to the tire performance. This requires identifying the key 
elements in the computational model of the tire and vehicle that are required to capture the major 
response characteristics. 


RELATED TASKS 

To accomplish the aforementioned technology needs, some related tasks need to be 
implemented before computational modeling of tires can have a significant impact on the tire 
design process. These are: 1) high-fidelity modeling of tire and tread pattern; 2) efficient 
computational strategies for rolling contact; 3) improved test methods, measuring techniques 
and scaling laws; and 4) validation of numerical simulations through a number of carefully 
designed benchmark tests. 
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Highlights of the Workshop 


Ahmed K. Noor 

Center for Computational Structures Technology 
University of Virginia 




OBJECTIVES AND FORMAT OF WORKSHOP 


The computational modeling of tire response represents one of the challenging applications of 
computational structural mechanics and has been the focus of intense efforts in recent years. Despite these 
efforts, major advances are needed in various areas of the technology before computational modeling can be 
routinely used in tire design. The objectives of the present workshop are (Fig. 1): to assess the state-of- 
technology in computational modeling of tires; and to identify current and future needs for further development 
of the technology. 

The workshop includes presentations and two panels. The presentations are included in the proceedings 
to illuminate some of the diverse issues, and to provide fresh ideas for future research and development. 


Objectives 


Format 

• Presentations 



Assess the state-of-technology in 
computational modelling of tires 

Identify future 
directions for 
research 


Panels 

- Panel 1 - Computational Needs for Tire Modeling and Design 
Moderators: R. Ridha and A. Noor 



- Panel 2 - Experimental Needs and Effective Coupling Between 
Experiments and Numerical Simulations. Moderators : J. Tanner 
and J. Luchini 


• Proceedings 


Figure 1 
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ASSESSMENT OF THE STATE-OF-TECHNOLOGY 


The first aspect of assessing the state-of-technology is to assess our understanding of the physical 
phenomena associated with tire response and durability. Some of the issues that affect these physical 
phenomena are listed in Fig. 2. These are: characterization of tire materials; footprint definition; noise, 
damage, wear, aging and failure mechanisms; friction modeling; accurate determination of operational 
conditions (viz., loads, strains, and temperature). 


• Understanding physical phenomena associated with 
tire response and durability 

- Footprint definition 

- Noise, damage, wear, aging and failure mechanisms 

- Friction modeling 

- Traction on both wet and dry surfaces 

• Characterization of tire materials 



Accurate determination 

- Loads and strains 

- Temperature 


of operational conditions 


Figure 2 
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ASSESSMENT OF THE STATE-OF-TECHNOLOGY (Cont'd.) 

Current Capabilities 


The second aspect of assessment of technology is that of current capabilities for computational modeling 
of tires (Fig. 3). These capabilities include computational material models ranging from micromechanical 
models for cord-rubber (flexible) composites and friction models to damage mechanisms, wear and failure 
models; modeling of the cross-sectional details using three-dimensional, quasi three-dimensional and two- 
dimensional finite element models; efficiency of currently-used computational strategies for handling 
frictional contact and rolling tire problem; and assessment of software systems currently used for modeling 
tire response and failure. 


Computational Material Models 

- Micromechanical modeling for cord-rubber 
(flexible) composites 

- Friction modeling 

- Damage mechanisms, wear and failure models 
Level of Details 

- Three-dimensional and quasi-three-dimensional 
models 

- Two-dimensional shell models 

- Other (detailed) models 
Efficiency of Computational Strategies 

- Frictional contact algorithms 

- Rolling tire problem 

Current Software Systems: MSC NASTRAN, 
ABAQUS, ANSYS, ADINA, . . . 




Figure 3 
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FUTURE DIRECTIONS FOR RESEARCH 


Three factors should be taken into account in identifying future directions for research (Fig. 4): 

1) Characteristics of future tires and their implications on design requirements; 

2) Future computing environment and computing paradigm; and 

3) Recent and future developments in other fields of comptuational technology which can be adapted 
to computational modeling of tires. 

Three of the important research tasks are: 

1 ) V alidation of numerical simulations and selection of benchmark tests for assessing new computational 
strategies and numerical algorithms. The standardized tests would provide a measure of confidence in added 
functional capabilities to existing codes, or in new codes; 

2) Treatment of uncertainties in material properties, geometry, loading and operational environment 
through probabilistic analysis, stochastic modeling and sensitivity analysis; and 

3) Understanding the relationship between the various response quantities and key performance 
measures for the tire (viz., durability and life, tread wear, noise, rolling resistance, vibrations, cornering, riding 
and handling qualities). 






Characteristics of future tires and their implications on 
design requirements. 

Validation of numerical simulations and effective 
coupling with experiments (benchmarks). 

Treatment of uncertainties in material properties, 
geometry, loading and operational environment 
(probabilistic analysis, stochastic modeling and 
sensitivity analysis). 




Impact of developments in other 
fields of computational technology 
(e.g., modeling of rigid composites, 
virtual manufacturing and IPPD). 


Impact of emerging and future computing 
environment (high-performance 
computers, multi-media workstations, 
advanced visualization technology). 


Figure 4 
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TOUCHDOWN DYNAMICS 


Samuel K. Clark 
Precision Measurement Co. 
Ann Arbor, MI 


INTRODUCTION 


Aircraft tire wear results from operating conditions are quite different from those encountered in land 
vehicles. One of the most important of these is touchdown, where the tire suddenly spins up from zero to a 
large angular velocity. This phenomenon is studied from both the analytical and experimental points of view. 
The analysis is basic, using elementary properties of the tire and wheel. It results in a new dimensionless 
description of the process. The experimental study consists primarily of small scale laboratory data, although 
limited full scale tire data is also presented. 

The data show increasing weight loss during touchdown as the dimensionless severity increases. 


'This work was conducted under USAF Contract F33657-92-C-2167. Mr. Paul Wagner was the USAF project engineer. 
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WHEEL AND AIRCRAFT PARAMETERS 


The analysis assumes that a tire and wheel are attached to the axle of an aircraft moving at forward 
velocity V x with constant sink rate V z , with the wheel at zero angular velocity to at time t=0 of initial contact. 
The tire of diameter D is attached to a wheel with brake rotor, all of which have polar moment of inertia I y . 
The vertical and horizontal contact forces are F z and F x = |iF z , respectively, where the friction coefficient p. 
is assumed constant during spin-up. 

The tire is assumed to be rigid torsionally, but to be a linear spring in the vertical direction, so that 
F z = K Z 8 Z where S z is vertical deflection and K z is the vertical spring rate. The contact area is controlled by 
tire deflection. 


v z 



Wheel and aircraft parameters 


F x = tire horizontal force 

V x = landing speed, assumed constant 

V z = sink rate, assumed constant 

l y = polar inertia of tire, wheel and brake 

K z = tire vertical spring rate 

ji = friction coefficient 

to = tire angular velocity 

D = tire diameter 

F z = tire vertical force 


Figure 1 - Wheel and tire parameters. 


10 



TOUCHDOWN ANALYSIS 


The analysis is based on evaluating the integral of contact area times sliding distance as a measure of the 
volume of tread material lost by abrasion and/or ablation. 

Assuming the volume of a tire is approximately proportional to the square of the diameter and the first 
power of the tire width, the resulting dimensionless loss of volume during touchdown is given in Fig. 2. This 
new parameter is called Touchdown Severity and is postulated to be a dimensionless measure of volume lost, 
as shown in Fig. 2. 

The friction coefficient p is not considered further in this analysis. 


Assumptions: 

(a) Vx> V z are constant during spin-up. 

(b) jli is constant during spin-up. 

(c) Contact patch area is proportional to tire deflection. 

(d) Wear is the integral of contact patch area x sliding distance. 


Result: 


volume lost _ 1 
volume of tire ” A 


; D 

> 

04 

W 

d 4 k z J 


V x = touchdown speed 
l y = polar moment of inertia 
D = tire diameter 
K z = tire vertical stiffness 


Touchdown Severity = 


V 2 I 

X y 


d 4 k 


A dimensionless number. 


Figure 2 - Touchdown Severity Index. 
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LABORATORY TESTING 


Based on the concept of dimensionless parameters, small scale tires were fabricated from Buna N (nitrile) 
O-rings mounted on nylon pulleys. These miniature wheels were mounted on an axle, which in turn could be 
forced into sudden touchdown to a rotating abrasive disc via a solenoid and a controlled timing circuit. By 
adjusting the point of touchdown on the disc from an inner position to the outer edge, touchdown velocity could 
be varied from 100 to slightly over 200 feet per second. By using inertia rings attached to the nylon wheel, 
polar moment of inertia could be varied. Using both of these as variables, the dimensionless severity index 
could be varied from values of approximately 1 0 up to 200. 

The timing circuit was adjusted to bring the wheel-tire into contact with the spinning abrasive disc for 
0.4 sec., a time sufficient to spin up, and then to retract it for 10 sec., allowing the wheel to spin down to rest 
by bearing friction effects. 

Ten touchdowns were used to make one complete set at a given condition. Weight of the tire-wheel was 
recorded both before and after the ten runs using a sensitive laboratory balance. 



Figure 3 - Small scale test apparatus for touchdown wear. 
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LABORATORY DATA ON BUNA N (NITRILE) 


Using the laboratory apparatus described on the previous page, data on Buna N (nitrile) rubber was 
obtained over a very wide range of values of the severity index, up a value of nearly 200. The values obtained 
were plotted as percent weight loss of the small tire due to ten touchdowns vs. severity index number. The 
data show a relatively smooth increase in weight loss as severity index increases. 
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Figure 4 - Percent weight loss vs. severity index. 




LABORATORY DATA ON NR AND NR-SBR COMPOUNDS 


Further data was obtained by fabricating small scale tires from compounds similar to those used in 
commercial aircraft tire treads. These were plotted against increasing severity index, and compared with the 
Buna N data previously obtained. As might be expected, the commercial compounds exhibited less wear. Two 
compounds were studied. The first, denoted No. 3 and plotted in Fig. 5 as circles, is a 100% NR compound. 
The second, denoted No.l3,isa65% NR-35 % S B R compound. They appear to be very close in their resistance 
to touchdown wear, but both are clearly more resistant than nitrile. 



Figure 5 - Weight loss at ten touchdowns vs. severity index. 
Three compounds. 
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ANALYSIS OF SHEAR STRESS 


It is known that at high values of the severity index cutting of the tread can occur. In order to examine 
the information available from this simplified analysis, the shear stress caused by touchdown is calculated by 
forming the ratio of the tangential force to the contact area. The resulting value of shear stress is probably as 
good a measure of cutting as is available. The resulting expression is primarily dependent on tire vertical 
stiffness K z , which in turn is controlled to a great extent by inflation pressure. This implies that cutting is most 
severe at high inflation pressures, but gives little other insight to the phenomenon, since there is in this analysis 
no simple means for defining the counter surface. 


Contact Area A = tcS^gjD = nV z tfwD 
Tangential Force F x = fiF z = fiV z t K z 

Shear Stress x x = F x /A = £ ~!=e 

xx n j wD 


Figure 6 - Shear stress as a measure of cutting. 
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SPIN DOWN TEST ANALYSIS 


In order to utilize existing spin down data, an analysis was conducted on the problem of a spinning tire, 
held at zero forward velocity, allowed to drop onto a rigid fixed surface. The result showed that the ratio of 
tread loss to tire volume is controlled by the same dimensionless severity index as obtained previously. In this 
particular case, a modified spin down experiment was considered in which the tire deflection was limited to 
preset value 8 by the use of air cylinders simulating aerodynamic lift. The specific value of 8 does not appear 
in the resulting volume loss. The rotational velocity decreases linearly with time in this analysis, as shown 
in Fig. 7. 


Modified Spin Down Test 


8 = allowable tire deflection 
Other assumptions same 


Yk.i 

v t -m 





Figure 7 - Volume loss and spin down speed. 
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SPIN DOWN TEST SPEED 


Full-scale spin down tests were conducted to evaluate tire wear under this condition. Instrumentation 
was installed to measure wheel and tire rotational speed during the test, and typical values of this data are shown 
in Fig. 8. The wheel speed decreases almost linearly with time, and this is in good agreement with the result 
of the analysis shown in Fig. 7. 


Slip vs. Time 



Figure 8 - Wheel speed vs. time on spin down test. 
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WEIGHT LOSS VS. TOUCHDOWN SEVERITY 


Using data from the full-scale drop tests, expressed as percent of tire weight for ten touchdowns, and from 
the small-scale NR and NR-SBR laboratory tires, a single combined plot of data can be produced as a function 
of touchdown severity. This shows that the resulting loss values are similar in magnitude, although the full- 
size data exhibits some scatter. The mass of the small tire was 2.4 grams, while that of the full-size tire was 
over 18,000 grams. Both bias and radial tires were used in the full-size test plan. 



TOUCHDOWN SEVERITY, V^I/D 4 K Z 

Figure 9 - Weight loss for both small-scale and full-size tires vs. touchdown severity index. 
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SEVERITY INDEX 


The touchdown severity index may be useful as a measure of the damage (abrasion, ablation or cutting) 
which is experienced by a tire during spin-up. Little reliable full-size tire data is available as yet to assess the 
validity of such predictions, but it is hoped that such data may be obtained in the future. 

This index may be interpreted as the ratio of the spin-up kinetic energy to the potential energy of the 
deflected tire, as shown in Fig. 10. 

i /(n 2 K)- K - E - of Tire & Wheel Assembly 
V [ u * z) = P.E. of Deflected Tire 


Figure 10 - Touchdown severity index. 
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REQUIREMENTS FOR ENERGY BASED CONSTITUTIVE 
MODELING IN TIRE MECHANICS 


John R. Luchini, Jim M. Peters and Will V. Mars 
Cooper Tire and Rubber Company 
Findlay, Ohio 


ABSTRACT 


The history, requirements, and theoretical basis of a new energy based constitutive model for 
(rubber) material elasticity, hysteresis, and failure are presented. Energy based elasticity is handled by 
many constitutive models, both in one dimension and in three dimensions. Conversion of mechanical 
energy to heat can be modeled with viscoelasticity or as structural hysteresis. 

We are seeking unification of elasticity, hysteresis, and failure mechanisms such as fatigue and 
wear. An energy state characterization for failure criteria of (rubber) materials may provide this 
unification and also help explain the interaction of temperature effects with failure mechanisms which 
are described as creation of growth of internal crack surface. 

Improved structural modeling of tires with FEM should result from such a unified constitutive 
theory. The theory will also guide experimental work and should enable better interpretation of the 
results of computational stress analyses. 
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INTRODUCTION 

Background 


This essay is part of our current research on "an energy based theory for three-dimensional 
elasticity, hysteresis, fatigue and wear." The word "essay" has its roots in the Latin verb for "to try." 
This presentation documents our current efforts, but not a final form for others to use. We hope that by 
presenting this admittedly preliminary discussion, we may stimulate other experts to aid in the search 
for an improved three-dimensional constitutive model for rubber. These others may offer insights as 
well as rigor to the concepts presented here. This three-dimensional constitutive model for rubber 
would ideally encompass three-dimensional elasticity as well as criteria for material damage or failure 
including hysteresis, fatigue and wear. 

The objective of this research is to develop a new, unified set of constitutive relations for an 
energy based material theory for the materials used in our product. Our emphasis on material property 
modeling was initially directed toward identifying better ways to provide input to the FEM computer 
programs and better ways to use material properties to interpret the output of the FEM programs. 

In an article [1] describing the fate of the ship 'Titanic" and the space shuttle "Challenger," a 
common theme was identified. In 1912, no one understood how the Titanic's steel would react when 
subjected to 31°F water temperature. The article claimed that many technological disasters could be 
traced to a common warning: know your materials. 

In the last decade researchers in tire mechanics have developed significant insights and 
understanding of the mathematical material model needed by the FEM technology [2-8]. We have 
achieved a level of technical sophistication with the FEM technology such that shortcomings in the 
laboratory measured material properties showed as a limitation in the utilization of FEM. 

Our current work on an energy based material theory started as a result of a meeting among 
personnel from research and materials development to discuss the differing requirements from testing 
for properties needed to use FEM, and to develop compounds. We saw that the current test results 
were distinct but that they should have some common basis which would allow a single set of tests to 
provide information which could be used for both FEM and materials development. 

We realized that a single underlying theory used in the creation of new material tests would not 
need separate tests for each application. The required test(s) may be more complex than any single 
current test but would provide a deeper understanding of the material being used or developed and 
would eventually lead to simpler and more distinguished tests. 


Thesis 

Energy theories for elasticity of a material are based entirely on the capacity of the material to 
generate, tolerate and contain energy as described in the Principle of Virtual Work and Castigliano's 
Theorem. The fatigue damage limit of a material might be entirely controlled by the capacity of the 
material to generate, tolerate and contain entropy. We propose a unified formulation of rubber behavior 
based on the Conversion of Energy Law and the Increase in Entropy Principle. 
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PROBLEM STATEMENT 
Generalized Unified Three-Dimensional Theory 
of Elasticity and Fatigue 


Energy based elasticity theories for the three-dimensional constitutive model for rubber should: 

- provide a theory to guide experiments and analysis of results 

- define inputs for FEM characterization of the material 

- provide material based criteria to use in FEM evaluations 

- combine what we know and what we need into a single framework 

- be based on scalar quantities such as the energy in the material. 

There are energy based elasticity theories, but we wanted to unify these with an energy based 
theory for hysteresis and accumulated damage. A unified energy theory should express theoretical 
relationships among material hysteresis, one-dimensional material modulus, three-dimensional elasticity 
modulii, and material fatigue. This theory must be applicable to rubber and should be consistent with 
SED fatigue theories, rubber incompressibility, and the DIOSS hysteresis theory. This model should 
include SED and temperature relationships to enhance the interpretation of FEM results. The 
constitutive model should be consistent with energy based elasticity as noted above and should also 
provide a theory for 

- an "energy capacity" of the material related to the energy input 

- the cycle of energy the material must survive 

- the hysteresis damage which may be accumulated before failure. 

For several years FEM studies have used and displayed a scalar quantity called "strain energy 
density" (SED) which "summarized" the stresses and made tire to tire comparisons simpler. While this 
provides an engineer with a simple graph to analyze, it does not interpret how different rubber 
compounds will react to similar levels of strain energy nor how many cycles a specific compound can 
tolerate at operating temperature. Thus, a better model for the material was needed but we were unable 
to find such a model in the technical literature. 
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PROBLEM STATEMENT 
Assumptions/Constraints for Elasticity 


The following is a list of the primary assumptions for the proposed material model: 

1. Uniaxial tension and compression stress-strain behavior which has been observed and 
modeled must be derivable from the general three-dimensional theory. Experiments performed in one 
direction must provide data for modeling the material behavior in other directions as well. 

2. Shear stress-strain behavior is isotropic. 

3. Shear stress-strain behavior is independent of the normal stress-strain behavior. This is 
based on data and analysis in a paper by Meineke and Porter [9]. 

4. Volume is constant, even at large strains. Large distortions must be accommodated by the 
theory, and initial volume and final distorted volume must be the same. Note that Poisson's ratio, 
which is defined for infinitesimal strains, will be a function of strain. Based on this fact, the classically 
defined Poisson's ratio should not be an input at all. 

5. Normal, volumetric, and shear elasticity parameters for the model may be used to derive 
equivalent values (at zero strain) for the classical elastic constants: Young's modulus, E, the shear 

modulus, G, and the Poisson's ratio, v. The derived values should correspond to the experimentally 
determined values in the literature. 

6. Normal, volumetric and shear elasticity in the model may be approximated as linear 
functions (single values) for low strains. However, for larger tensile, compressive, or shear strains, 
they must be nonlinear functions. 

7. The elasticity model must provide for or permit functional relationships of stresses and 
strains with temperature. Frequency effects, or other parameters, which are determined to be important 
must be handled in a way that is consistent with the hysteresis and failure theories which are under 
parallel development. 
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PROBLEM STATEMENT 

Assumptions/Constraints for Hysteresis and Fatigue 


A paper presented at the Tire Society last year pointed out that there is no uniquely good way to 
analyze and interpret tire three-dimensional stress analysis problems [2]. To develop a new fatigue 
theory, we are considering available theories of polymer fatigue from a distortional energy basis. This 
has led us in the direction of an entropy capacity fatigue model. The theoretical model we desire will: 

1 . View hysteresis as a result of an internal (dry) friction process [3] which creates an 
accumulation of entropy within the material micro-structure. 

2. Blend with, or extend, the DIOSS theory for hysteresis related to distortional energy to 
include fatigue and elasticity. This can be modeled as: H=B*exp(-rT)*(A0SS A 2)/(l+P0SS). 

3. Demonstrate that higher hysteresis within the material reduces crack growth rate due to 
energy dissipation. 

4. Have material properties or parameters which can be used in tire stress analysis. There are 
theories which may work but are incompatible with the other assumptions used in tire analysis. 

5. Be implemented into FEM and coupled with stress/strain analysis for durability prediction. 

6. Avail itself of straight-forward laboratory tests which can be used directly for materials 
development, and provide data for material properties or parameters. 
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PROBLEM STATEMENT 
Reasonable Experiments to Obtain Parameters 


The motivation for needing a new constitutive relationship is the unique behavior of filled, 
vulcanized polymers. While the points listed in the previous page must be part of the foundation for the 
new three-dimensional model, they are not sufficient to specify die experiments which define the 
parameters for the material properties. This parametrization will probably require unique tests and 
measurements to define material parameters. 

Energy storage modes within the material would include: strain energy; thermal energy; 
chemical energy stored in bonds; and a "surface energy" which increases as surface is 
formed as in crack growth or fracture. 

Entropy capacity of the material would include factors such as: the ability of the chemical bonds 
to stretch, bend, or rotate; volume of the material (for "entropy density"); and the spatial 
distribution of material components of different entropy capacities including both structural 
features and dispersion. 

Characterization of rubber compound behavior in the past has utilized various tests. Future 
experimental analysis techniques must be consistent with these results. In addition, rubber compounds 
used in tires exhibit characteristics which the new constitutive model must include. The features of the 
model must: 

- provide for a nonlinear relation between stress and strain 

- provide for anisotropic behavior, and describe stress and strain uniquely in different directions 

- provide for incompressibility and permit distortion with no volume change. There is a 
temperature volume change effect. 

- possess a firm theoretical foundation, be experimentally verifiable and compatible with 
existing FEM theory. 
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ENERGY BASED ELASTICITY THEORY 
Elasticity From Energy 


The Mooney-Rivlin material model quantifies rubber modulus in terms of two tensor invariants. 
One of these relates to the energy stored in the material by compression and the other energy stored by 
shear. It is generally agreed that all energy stored by a compressive stress on an incompressible 
material is elastic, that is, fully recoverable. However, some of the energy input as shear will not be 
recovered. Experimental measurements of the Mooney-Rivlin constants for carbon black filled rubber 
show different values in extension, compression and shear. In other words, the "constants" depend on 
the strain conditions and would require complex computer coding to obtain a result which is both 
physically meaningful and mathematically stable. 

The use of a scalar quantity such as SED or energy is a common approach in constitutive 
material modeling because it permits the development of an unambiguous three-dimensional theory 
from a theory developed in one dimension. We have noted that the existence of constant values for E, 

G and v imply a volume change for the material under stress. We have conducted laboratory 
measurements to confirm the non- tire industry research published in 1993 [10] on a new constitutive 
model for rubber using a variable for Poisson's ratio. The assumption that the shear modulus and 

material volume must remain constant implies that the modulus, "E", and/or Poisson's ratio, "v", must 
generate a nonlinear relation between stress and strain. 

Nonlinear modulus theories which are not based on scalar quantities can develop mathematical 
difficulties or predict nonphysical results when extended into three dimensions. One such case is the 
use of nonlinear rubber modulus in the current Cooper FEM program. This created directional 
properties in the elements under stress which were not seen in real tires and which caused 
mathematically unstable computations. The incompressibility material model in the current Cooper 
FEM program must be computationally complex in order to prevent these mathematical difficulties. 

The "Directional Incremental Octahedral Shear Strain" theory ('DIOSS") quantifies rubber 
hysteresis in laboratory experiments and provides parameters for use in FEM analyses. This theory 
uses a scalar quantity. Octahedral Shear Strain (OSS), to compute the quantity of energy lost to 
hysteresis. The research completed shows that the OSS is related to the distortional energy 
mechanically put into the material, some of which is lost to hysteresis. 

A combination of the Mooney-Rivlin tensor invariant related to the energy stored in the material 
by compression and the DIOSS scalar quantity for energy stored by shear and lost to hysteresis could 
form the foundation for a new three-dimensional constitutive model. The model would retain the 
hyseretic response from our DIOSS theory. The new constitutive model could also have terms for the 
elastic response which would replace current elastic modulus in FEM. 

The features we desire in the theory should be based on tensor invariants in order to provide a 
nondirectional material model. This is desirable in order to be both physically meaningful based on the 
observed experimental data as well as providing for a mathematically stable solution. 
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ENERGY BASED ELASTICITY THEORY 
Thermo-Mechanical Interaction 


We have noted that a mechanical cycle puts energy into the material and then removes it. While 
we can remove enough energy to return the material to its original energy density, some entropy 
increase will persist in the material after the cycle is complete. 

In several papers [2-8] the thermo-mechanical interactions which occur in tires have been 
modeled with sequential applications of FEM stress analysis, a stress-strain to hysteresis thermal 
energy estimate, and a thermal analysis. Some of these approaches allow for iterative solutions, but 
none directly allow for interaction of the thermal and mechanical energies. A diagram, shown below, 
describes a way these energies may interact. 
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ENTROPY BASED DAMAGE THEORY 
Damage Criteria Aids FEM Interpretation 


The general answer to the need for a three-dimensional elasticity theory should be in terms of 
energy. Elastic strain energy is one form of energy and hysteresis is another form. There are problems 
with this approach in terms of the mechanical analysis available from FEM. We are able to compute the 
response of the material in three-dimensional space, but the computation of energy (either strain energy 
or hysteretic energy) eliminates the directional information. We need a computation which can be done 
within FEM before the directional information is lost. 

An analog of the computation required is the "wind-chill index" from the weather reports. 
Temperature is a scalar quantity similar to energy. However, your body is actually responding to heat 
transfer (energy flux) when you say you are cold. A correspondence between a person's level of 
comfort in a given "still air" temperature and the current temperature and wind velocity, related via heat 
transfer coefficients, is the "wind-chill index." 

We need to identify what the material is actually responding to at a given level of energy. Then 
we may be able to create one or more indices which can be computed within the FEM program to 
correspond to the material's potential discomfort. Since we are dealing with tensors rather than 
vectors, we expect that at least two indices would be needed. 

As an example, consider an "elastic strain-energy index" and an "hysteretic energy index" 
which could have critical values for various failure modes. Another set of indices could be a dilational 
energy index" and a "distortional energy index" which could have critical values for energy storage and 
energy dissipation, respectively. Unfortunately, we do not know which set of indices are desirable to 
compute, nor how to compute these desirable indices in a way which accurately represents the 
material's level of discomfort. 
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ENTROPY BASED DAMAGE THEORY 
Thermal Hysteresis and Entropy 


The concept which has emerged has variously been called a "material wind-chill index," the 
"material energy limit," the "energy capacity," or the "hysteretic energy index." The hysteretic energy 
index would be a material property which describes the ability of a material to tolerate hysteresis. A 
term we have coined, the "thermal hysteresis" would be the heat energy left over after a thermal loading 
has been applied to a system and then removed. We now think the idea we have been describing is 
actually entropy generation within a material. 

For example, consider a slab of rubber heated from 75°F to 150°F and then cooled once again 
to 75°F. The difference between the thermal energy input to the rubber sample and the thermal energy 
returned after cooling would be zero for an ideal material. However, the heating process is not 
completely reversible and the difference between the thermal energy input and the thermal energy 
returned would be the thermal hysteresis of the material. Based on the first law of thermodynamics the 
above may be written as: 

Thermal energy input = thermal energy output + thermal hysteresis 

This thermal hysteresis process is the same as the principle of increase in entropy from the 
second law of thermodynamics. While entropy may be conserved during idealized reversible 
processes, entropy is generally a nonconserved property. There is no such thing as the conservation of 
entropy principle, and it increases during all actual processes. 

Heat transfer to a system increases system entropy along with the energy increase. Heat 
transfer from a system decreases both the energy and the entropy of the system. Entropy change as a 
result of reversible heat transfer is called entropy flow and does not involve any entropy generation. 

Irreversibilities such as friction, fast expansion or compression, and heat transfer through a 
finite temperature difference always cause the entropy to increase. It is hypothesized that the material 
damage is a result of the entropy increase and/or can be quantified as proportional to this thermal 
entropy increase. 

Energy theories for elasticity of a material are based entirely on the capacity of the material to 
generate, tolerate, and contain energy, as described in the principle of virtual work and Castigliano's 
Theorem. The fatigue damage limit of a material might be entirely controlled by the capacity of the 
material to generate, tolerate, and contain entropy rather than energy. 

Thermal energy input by itself is not assumed to cause material damage, although it can change 
the physical properties of the material and thus affect the loading conditions. The damaging portion of 
the thermal energy, the "thermal hysteresis," is the entropy increase which results from chemical 
changes in the material. This entropy increase would be a function of time, and could also be used to 
explain the effects of aging on rubber. 

We hypothesize that, for a given application, there is a "material limit" for entropy generation, 
instead of energy content. In such a material constitutive model, the material's life is determined by the 
quantity of entropy generated throughout its history. When this quantity exceeds the limit, it is 
considered to have failed. This principle should be reconcilable with accumulated damage theories such 
as Miner's Rule. 
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ENTROPY BASED DAMAGE THEORY 
Thermal Hysteresis and Entropy (Cont'd.) 


If we could quantify a property of the rubber for the entropy generation rate as a function of 
strain and temperature, we should be able to integrate over time to find the total entropy generated. 
This entropy generated and accumulated in the material would quantify the total damage done in the 
cycle of strain and temperature. By these hypotheses, the energy content at a given time and the total 
accumulated entropy would determine the elasticity and life of the material. The damage done to the 
material would be equivalent to that done by heat and could be called the thermal hysteresis or 
accumulated entropy generation. 
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ENTROPY BASED DAMAGE THEORY 
Hysteresis and Fatigue Damage 


This paper proposes that fatigue damage is the accumulation of entropy. The discussion here 
also implies that thermal hysteresis by itself could be used to predict durability. The "energy limit" of a 
material we initially described is now seen as the amount of thermal energy needed to increase the 
entropy beyond the material's tolerance and thus cause material failure. 

This "energy/entropy limit" is still an imprecise concept and we have no straight-forward 
laboratory process to measure it. However, there is a strong analogy between what happens to the 
rubber mechanically and what happens thermodynamically during a fatigue process. The following 

diagram depicts the analogy between heat energy, Q, and mechanical stress, o. 
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CONCLUSION 

Hysteresis, Thermal Energy and Entropy 


While there are good energy based elasticity theories, we have been unable to find or develop a 
good entropy/energy based damage theory. We, therefore, are proposing "HEAT," the hysteretic 
entropy accumulation theory. HEAT claims that entropy generation occurs due to irreversible 
microscopic structural changes, at the elastomer network level. Energy loss occurs due to irreversible 
macroscopic structural changes, at the filler level. While the friction generated macroscopically would 
not damage the rubber, irreversible changes in the elastomer network would damage the rubber. 
Energy lost due to macroscopic irreversibilities shows up as heat, part of which ends up causing 
microscopic irreversibilities. 

This model would encompass chemical interactions that occur in the rubber. Chemical 
interactions (changes to the molecular network) would be the source of entropy generation, just as 
friction between filler particles is the source of energy loss. 

One reason that hysteresis is historically related to fatigue damage may be due to the heat 
generation which results from mechanical losses. In a "real" material mechanical hysteresis presents 
itself as a temperature rise. This temperature rise is a thermal energy input, and as that thermal energy 
input dissipates into the environment the material would accumulate an entropy increase. It is our 
hypothesis that the entropy is either the cause of the damage or the result of the damage, but in either 
case it provides a measure of the damage. The material property which we seek is one which can be 
used to quantify the material's tolerance of the entropy accumulation. 
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CONCLUSION 
Model Form and Flow Chart 


A schematic idea of how the HEAT fatigue/damage accumulation process would operate is 
shown below. This model not only addresses the mechanisms of fatigue, but could also be used to 
describe rubber aging. 
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A KINEMATICALLY DRIVEN ANISOTROPIC VISCOELASTIC CONSTITUTIVE 

MODEL APPLIED TO TIRES 


Arthur R. Johnsont, John A. Tanner^, and Angela J. Mason* 
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INTRODUCTION 


Aircraft tires are composite structures manufactured with viscoelastic materials such as carbon black 
filled rubber and nylon cords. When loaded they experience large deflections and moderately large strains. 
Detailed structural models of tires require the use of either nonlinear shell or nonlinear three dimensional 
solid finite elements. Computational predictions of the dynamic response of tires must consider the 
composite viscoelastic material behavior in a realistic fashion. In this extended abstract we describe a 
modification to a nonlinear anisotropic shell finite element so it can be used to model viscoelastic stresses 
during general deformations. The model is developed by introducing internal variables of the type used to 
model elastic strain energy. The internal variables are strains, curvatures and transverse shear angles 
which are in a one-to-one correspondence with the generalized coordinates used to model the elastic strain 
energy for nonlinear response. A difference-relaxation equation is used to relate changes in the observable 
strain field to changes in the internal strain field. The internal stress state is introduced into the equilibrium 
equations by converting it to nodal loads associated with the element's displacement degrees of freedom. 

In this form the tangent matrix in the Newton-Raphson solution algorithm is not modified from its form 
for the nonlinear statics problem. Only the gradient vector is modified and the modification is not 
computationally costly. The existing finite element model for the Space Shuttle nose gear tire is used to 
provide examples of the algorithm. In the first example, the tire's rim is displaced at a constant rate up to a 
fixed value. In the second example, the tire's rim is enforced to follow a saw tooth load and unload curve 
to generate hysteresis loops. 


t MS240, Army Research Laboratory, Computational Structures Branch. 
$ MS497, Structural Dynamics Branch. 
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SPACE SHUTTLE NOSE GEAR TIRE LOADING AND UNLOADING DATA 


Loading and unloading a Space Shuttle nose gear tire results in a hysteretic load displacement curve. 
The curve shown below was taken from Ref. 3. It is clearly seen that a significant amount of energy is not 
recovered during unloading. The objective of this effort is to develop a computational method to model 
general time dependent deformations of composites using stress and strain measures similar to those used 
to model the elastic (time independent) material behavior, that is, to make maximum use out of the existing 
form of the nonlinear finite element code and the energy method which it is based upon. The model can 
then be used to predict the hysteretic behavior of tires so as to improve their design for enhanced 
performance and life. 


The work presented in this extended abstract is limited to a brief description of the computational 
method developed and of its use to simulate tire loading and unloading. 


Vertical 

load, 

kips 



Deflection, w, in. 


Load deflection data for Space Shuttle nose gear tire. 
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NONLINEAR SHELL ELEMENT AND TIRE MODELING CODE 


A shell element capable of modeling geometrically nonlinear deformations of laminated composites 
developed by Noor, et al 1*2 [ s shown below. It is the element used in NASA's tire footprint finite 
element code. This shell element has nine displacement nodes with five variables at each node and four 
stress nodes with eight variables at each. 


The NASA tire footprint code^ allows the user to generate a tire stress analysis model based on shell 
theory. The code can then simulate the tire inflation and enforced rim displacements with the tire pressing 
against a frictionless flat surface (for a surface with friction see Ref. 4). 



x Stress Node 

( T| £ ) Nondimensional Element 
Coordinates 


43 



ELASTIC CONSTITUTIVE EQUATION 

The stress strain relation used to model the composite shell is given below. The shell's inplane 
stresses are N s , Ng and N S q , moments are M s , Mg and M S g ; and transverse shear forces are 

Q s and Qg . The shell's inplane strains are £ s , £q and 2 £ s q‘, curvatures and twist are 
k s , Kq and k s q ; and transverse shear strains are 2 e s 3 and 2 Eq^- 



We abbreviate the constitutive matrix as follows: 


C f 0 £ 

f T d 0 - k 

0 0 Cy 7 

where, 

£ 

K ► = the Sanders-Budiansky nonlinear strains^’^. 

y. 
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PLATE ELEMENT STRAINS AS INTERNAL VARIABLES 


Johnson, et al^-14 recently demonstrated that the three - dimensional strain dependent viscoelastic 
stresses measured in rubberlike materials can be modeled using either time dependent reference shapes (an 
internal solid theory) or an auxiliary set of time dependent strain variables (a kinematic variable theory.) 
Here, we follow the second method mentioned above. A set of internal strain variables and a constitutive 
relation are introduced to specify the magnitude of the viscous stresses. The continuous time dependent 
deformation of the solid is approximated by discrete strain steps in time. A step strain relaxation finite 
difference equation then drives the internal strain variables magnitudes at each discrete step of the physical 
system and also partially relaxes the internal strain variables between the discrete steps. 


Under the assumption of time - strain separability the viscoelastic stresses in the shell element at any 
time can be written as follows: 
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where N ( t ) , M ( t ) and Q ( t ) are the viscous components of N(t), M(t) and Q(t) 
respectively. It is assumed that the constitutive model for the viscous stresses is given by: 
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where e (t), k (t) and y ( t ) are internal kinematic variables (i.e., they are strain measures not 
stress measures.) The components in the constitutive matrix are assumed to be constants. 
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KINEMATIC EXCITATION OF THE INTERNAL VARIABLES 


During a general time dependent deformation of a solid the stresses measured in any given strain 
state will be different than the long term (relaxed) stresses in that strain state. Molecular theories of 
viscoelasticity (see Ref. 15) use a measure of the solid's deformation per unit time to drive the internal 
variables during an increment in strain. The internal variables such as molecular entropy and molecule 
length are then relaxed during the time between strain steps. This technique (see also Ref. 12) is applied to 
the shell's internal variables as follows: 
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Deformed Solid's Strain 


o o 


-o 



-e 

t Time 

3 



The finite difference update equation used to drive and define the internal variables is 


e j + 1 ~ 




e~ At / T 


The relaxation factor, T, is selected by fitting the constitutive model to step-strain relaxation data. This is 
an explicit update and results in an implementation with minimal changes to the finite element code. It 
does, however, require the use of small time steps (and several executions at different time step sizes) to 
maintain accuracy. 
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EQUILIBRIUM EQUATIONS 


In this model the viscous stresses due to the internal (state) variables at any given time are converted 
to applied forces at the displacement nodes. The finite element equilibrium equations are then modified by 
changing the load vector. The total stress at any time within the shell element is the sum of the elastic and 
viscous stresses. Equilibrium is then achieved between the stresses and the loads due to applied pressure, 
supports, and contact. 


Using the notation of Noor, et al*' 4 and the equations introduced above, one has at any time, t, the 
viscous stresses at the stress nodes given by: 

KHKFJ 

f * | at the displacement nodes using the interpolation functions. 


These stresses are converted to forces 
Again, using the notation of Noor, et al the following expression results: 

iT 


KHFK} 


where [ j = the linear strain displacement matrix. The nonlinear force contributions due to the viscous 
stresses undergoing large displacements is not computed in this effort. 


The Newton-Raphson equations used to find equilibrium are modified by subtracting j f* j from 
the gradient vector determined in the Taylor expansion of the equilibrium equations. It becomes: 


s £ + m n£ 
0 


s £ + m n£ 

0 


'Ah' 


rh 

0 


« 

Ax 

> = < 

* 

rx — f • 

^ c 

£c 
e . 


AX 


rX 


where f ^ = the linear flexibility matrix, m n ^ = the nonlinear strain displacement matrix, q c = a matrix 
associated with contact, r c = a matrix for regularization of the functional for contact, e = a penalty 

parameter. Ah , Ax , AX = stress resultants, nodal displacements, and Lagrange (contact) multipliers, 
rh , rx , rX = gradient expressions from the Taylor expansion of the equilibrium equations. 
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SIMULATION OF RELAXATION AND HYSTERESIS 


The NASA tire footprint finite element code was modified to implement the model described abovet. 
Since experimental data was not available to calibrate the viscous material law a number of trial simulations 

were made with different relaxation times, t, and with different constitutive matrices, c* J, to determine 

sample values for use in a computational simulation. We used t = 1 . 0 sec for the relaxation of the 

internal variables, and [ c* ] = [ c ] . The viscous moduli were assumed equal to the elastic moduli since 

the internal variable update relation is explicit in time. Recall: 


"j + 


i = [<$ + ( 


- 1 " £ J 


-4t / T 


A number of trial computer executions were used to determine an appropriate time step size. That is, a 
numerically determined time step, small enough for convergence, was used; namely 
At = 0. 0625 sec. 


In the first simulation the rim displacement was ramped to 1.2 in. over a time of 5.0 sec. The 
displacement was held constant. The computed force-time curve is shown below. Both the total stress 
(viscous) and the elastic (static) stress are shown. 



Time (sec) 



t The authors wish to thank Ms. Jeanne M. Peters of the Center for Computational Structures 
Technology, University of Virginia, NASA Langley Research Center for her assistance in making the 
necessary changes to the finite element code. 
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SIMULATION OF RELAXATION AND HYSTERESIS (Continued) 


The second simulation consisted of forcing the tire rim from 0.0 in. to 1.2 in. over a period of 5.0 
sec and then back to 0.0 in. during the next 5.0 sec. This simulation produced the hysteretic loop shown 
below. The load displacement curve for the static loading (without any viscoelasticity) is also shown. 


G 

•H 


1.5 


4 -> 

G 

0 ) 

e 

Q) 

U 

(0 

m 

■H 

Q 

e 

•H 

PS 


1 


0.5 



Time (sec) 


u 


Rim Displacement (in) 


In the third simulation the same rim displacement rate was used, but the rim was displaced to 1.8 in. 
before it was unloaded. It was then displaced to 0.9 in. The total (viscous) and elastic (static) stresses are 
shown below. 


Load - Unload - Load 
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SIMULATION OF RELAXATION AND HYSTERESIS (Continued) 


A detailed look at the Contact Force vs Rim Displacement graph shown below indicates that the 
contact force went to zero when the rim was "unloading". This is due to the fact that the "tire" did not 
have time to fully recover during the unloading. It separated from the contact surface. As is shown 
below, while the contact surface continued to drop and rise the tire recovered slightly and contact was 
made at a smaller rim displacement upon reloading. 



0 0.1 0.2 0.3 0.4 0.5 


Rim Displacement (Ins) 


SUMMARY 


A method was presented in which the mixed nonlinear shell element of Noor, et all >2 was use d to 
simulate viscous deformations of a Space Shuttle nose gear tire during a vertical loading and unloading 
against a frictionless surface. The viscous constitutive model introduced is anisotropic and is defined by 
internal variables in one-to-one correspondence with the generalized coordinates used in the shell's 
constitutive model. The computational simulations were limited to a demonstration of the algorithm since 
experimental data was not available to calibrate the viscous model. The algorithm appears promising in 
that it can be modified to include more complicated time dependent material behavior than was 
demonstrated above. That is, additional internal variables can be introduced if several relaxation times are 
needed to model the material's behavior. 


50 




REFERENCES 


1. Noor, A. K. and Hartley, Nonlinear Shell Analysis via Mixed Isoparametric Elements, Computers & 
Structures, Vol. 7, pp 615-626, 1977. 

2. Noor, A. K., Anderson, C. M., and Tanner, J. A., Exploiting Symmetries in the Modeling and 
Analysis of Tires, NASA Technical Paper 2649, March 1987. 

3. Kim, K. O., Tanner, J. A., Noor, A. K., and Robertson, M. P., Computational Methods for 
Frictionless Contact With Application to Space Shuttle Orbiter Nose-Gear Tires, NASA Technical Paper 
3073, May 1991. 

4. Tanner, J. A., Computational Methods for Frictional Contact with Application to the Space Shuttle 
Orbiter Nose-Gear Tire, PhD Dissertation, George Washington University, September 1992. 

5. Sanders, J. L., Nonlinear Theories for Thin Shells, Q. Appl. Math., vol. 21, no. 1, pp 21-36, 1963. 

6. Budiansky, B., Notes on Nonlinear Shell Theory, J. Appl. Mech., vol. 35, no. 2, pp 393-401, 1968. 

7. Johnson, A. R., Quigley, C. J., Cavallaro, C. and Weight, K. D., A large deformation viscoelastic 
finite element model for elastomers, in The Mathematics of Finite Elements and A pplications VII. Edited 
by J. R. Whiteman, Academic Press Limited, ISBN 0-12-747257-6, 1991. 

8. Johnson, A. R., Quigley, C. J., Weight, K. D., Cavallaro, C. and Cox, D. L., Inflation and deflation 
of a thick walled viscohyperelastic sphere, in Transactions of the Eighth Army Conference on Applied 
Mathematics and Computing, U. S. Army Research Office Report No. 91-1, 847-857(1991). 

9. Johnson, A. R. and Quigley, C. J., A viscohyperelastic Maxwell model for rubber viscoelasticity. 
Rubber Chemistry and Technology, 65, 137-153(1992). 

10. Johnson, A. R., Quigley, C. J. and Freese, C. E., Viscohyperelasticity, in Transactions of the Tenth 
Army Conference on Applied Mathematics and Computing, U. S. Army Research Office Report No. 93- 
1,235-256(1993). 

11. Johnson, A. R., Quigley, C. J., Young, D. G. and Danik, J. A., Viscohyperelastic modeling of 
rubber vulcanizates, Tire Science and Technology, TSTCA, 21, No. 3, 179-199(1993). 

12. Johnson, A. R. and Stacer, R. G., Rubber viscoelasticity using the physically constrained system's 
stretches as internal variables, Rubber Chemistry and Technology, 66(4), 567-577(1993). 

13. Johnson, A. R., Quigley, C. J. and Mead, J. L., Large Strain Viscoelastic Constitutive Models for 
Rubber. Part I: Formulations, accepted for publication. Rubber Chemistry and Technology, 67, (1995). 

14. Quigley, C. J., Mead, J. L. and Johnson, A. R., Large Strain Viscoelastic Constitutive Models for 
Rubber. Part II: Determination of Material Constants, accepted for publication. Rubber Chemistry and 
Technology, (1995). 

15. Doi, M. and Edwards, S. F., The Theory of Polymer Dynamics, Oxford, New York, 1986. 


51 




The Mode I Mechanical Crack Tip Stress Field in 
Hyperelastic and Incompressible Materials 


Claudia J. Quigley 
U.S. Army Research Laboratory 
Materials Directorate 
Watertown, MA 


53 




THE MODE I MECHANICAL CRACK TIP STRESS FIELD IN HYPERELASTIC AND 

INCOMPRESSIBLE MATERIALS 


Claudia J. Quigley 
U.S. Army Research Laboratory 
Materials Directorate 
Watertown, MA 


INTRODUCTION 


The finite deformation field of a plane strain Mode I crack in a hyperelastic and incompress- 
ible material was examined under the assumptions of small scale nonlinearity. Finite element 
analyses were performed for two different material laws, a Neo-Hookean material and a third 
order invariant of a Rivlin material. The numerical results for both materials were compared to the 
appropriate theoretical asymptotic solution. A local cavitation locus surrounding the crack tip was 
identified for the Neo-Hookean material. For the third order invariant Rivlin material, maximum 
values of the dominant stress component were found close to the surface of the crack, above and 
below the deformed crack tip. 
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PROBLEM DESCRIPTION 


Consider a thick infinite plate of hyperelastic and incompressible material, as shown in Fig- 
ure 1, containing a stationary crack. The crack faces were traction free and a uniform uniaxial 
stress, x, was applied at infinity. No body forces or thermal forces were present The stresses and 
deformation field in a small circular region encompassing the crack tip were examined to deter- 
mine the size and the shape of the nonlinear crack tip region. Crack tip field quantities were com- 
pared to the theoretical asymptotic solution for this problem. 
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Figure 1. A crack in a thick, infinite plate of hyperelastic and incompressible material. 
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APPROACH: SMALL SCALE NONLINEARITY 


The approach taken for this problem was that of small scale nonlinearity 1 , which is analo- 
gous to small scale yielding. Referring to Figure 2, consider our circular region of radius, R. The 
crack length within this region is also of length, R. Under small scale nonlinearity, the applied 
load is sufficiently small to confine the crack tip field to a small region surrounding the crack tip, 
compared to the crack length and the overall geometry of the body. Three separate regions of 
material behavior can be identified. Far from the crack tip, small geometry changes are present, 
and material behavior approaches linear elasticity. Here, the linear elastic Mode I asymptotic solu- 
tion was present The linear elastic region encloses a smaller, second region, where material 
behavior transitions from linear elastic to nonlinear elastic. This second region encompasses a 
third, smaller region, where finite deformations and large deformations are present, indicative of 
nonlinear elastic material behavior. The nonlinear asymptotic solution is dominant in this region. 





Nonlinear Region 
Transition Region 


Linear Elastic Region 


Figure 2. Small scale nonlinearity. 
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ASSUMPTIONS OF SMALL SCALE NONLINEARITY 


Under small scale nonlinearity, the J-integral 2 is assumed to be path independent around any 
suitable contour, T, which surrounds the crack tip, as shown in Figure 3. The J-integral can be 
expressed as 


J = 


\(W ni - Si u iA )dS 

r 


where W represents the strain energy, nj is the unit normal in the xj direction, Sj is the nominal 
traction vector, and Uj is the displacement vector. Here, dS is an element of arc length along T. 
The J-integral is also assumed to be equivalent to the energy release rate, so that the value of the J- 
integral is a direct function of the applied stress intensity factor, Kj. 


J = 


.2 U-'u) 

I 2\i 


Here, |i is the shear modulus, and the Poisson ratio, o, equals 1/2 due to incompressibility. 



Figure 3. Evaluation of the J-integral around a suitable contour. 
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BOUNDARY CONDITIONS 


Due to reflective symmetry about the xj axis, only the top half of the circular plate was mod- 
elled. Consistent with small scale nonlinearity, the Mode I linear elastic solution was applied as 
consistent nodal point forces along the outer radius, R, see Figure 4. The linear elastic solution 3 is 
given as 

K I 6 . . 0 . 30. , wo. 

a ll = ^ “ sm ^ sm — ) +o (r 1/z ) 


0 


. 0 . 30, 


1 / 2 , 


o 22 = J== cos^ (l + sin^siny) +o(r 1 "') 


0 30 . 0 


. 1 / 2 , 


°12 = J^ C0S 2 C0S T Sm 2 +0(r } 


in terms of the undeformed polar coordinate system (r,0).At sufficiently small loads, the linear 
elastic solution is found throughout the entire crack tip region. However, as the magnitude of Kj is 
increased, the nonlinear crack tip field began to evolve at the crack tip and grew in size and extent 
It is the stresses and deformations in the nonlinear crack tip region that will be examined next 



Mode I Linear Elastic 
Solution Applied Along 
the Outer Radius 


Figure 4. Boundary conditions applied along the outer radius, R. 
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FINITE ELEMENT MODEL 


The finite element mesh 4 was designed to sustain very large deformations at the crack tip, so 
that the nonlinear region could sufficiently expand to a size at which the nonlinear crack tip field 
quantities could be compared to the theoretical asymptotic solution. At these large deformations, 
numerical mapping methods begin to break down, and the elements begin to evert, making die 
nonlinear crack tip field difficult to resolve. Consequently, the finite element mesh. Figure 5, was 
composed of two regions, a coarse mesh in the immediate vicinity of the crack dp, encompassed 
by a refined mesh, where accurate field quantities could be found. 



-R Crack Face ^ R 

Crack Tip 


Figure 5. The finite element model. The coarse mesh in (a) is enclosed by the refined mesh in (b). 
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THE NONLINEAR ASYMPTOTIC SOLUTION FOR THE 
NEO-HOOKE AN MATERIAL 


For the first finite element analysis, a Neo-Hookean constitutive law was chosen. For this 
material law, the strain energy is expressed as: 

W = £(/- 3) 

where the principal strain invariant, I, is a function of the stretch ratios, Xj, 

/ = X* + X| + X 2 

Under the assumptions of small scale nonlinearity, the Mode I nonlinear asymptotic solution 5 for 
this material is: 


•01 3/2 x- 6 2 30. 

yj-ftrsin j + -f /i (lcos J - jCOSy) 


H 2 — 1 
X 22~4 fl r 


tjj-^^ci-cose) 


1/2 6 , _ 3/2 / .„• 30 b A 


e. 


y 2 ~ ar ' sin^ + r D/ * (dsm— - — srn^) 


x n~\ abr 1/2 S “I 


where a, b, and d are load parameters. The defamation field is given by y lt and the Cauchy 
stresses are represented as ty. Substituting this asymptotic solution into die definition of the J- 
integral permits a direct solution for the load amplitude parameter, a, in terms of Kj. 


a = 


*ijjt 


From the finite element analysis, the remaining two load parameters were found. 


b — 1.6 


when Jappiied = 0.065 MPa-m 


d~ 
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RESULTS FOR THE NEO-HOOKEAN MATERIAL 


The graph below is a normalized plot of T 22 48 a function of undeformed radial distance 
from the crack tip. These results were taken at an angle of 71.25 degrees. Far from the crack tip, at 
large values of r, the linear elastic asymptotic solution is found. As the crack tip was approached, 
the finite element results transition to the nonlinear elastic solution for this problem. Figure 6 is 
only one example in which the presence of the nonlinear elastic solution was verified. To see addi- 
tional results, the reader is referred to reference 4. 




= o.056 
[LR 



^ applied^ V" 


Figure 6. Logarithmic plot of T 22 as a function of radial distance at 0 = 71.25 degrees. 
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FAILURE MECHANISM 


A potential failure mechanism in this material is cavitation of a pre-existing microvoid. For 
a Mode I crack, representative of a nonhydrostatic load condition, a cavitation criteria*’ has been 
expressed in terms of the principal Cauchy stresses, T;, as 

(4t 3 ~ T l “ (4t 2 “ X 1 " V ^ (4t l ~ x 2~ V “ 125 ^ 3 = 0 

This criterion can be applied to the mechanical crack tip stress field to determine a surface sur- 
rounding the crack tip where cavitation is likely to occur. This region is outlined by the dashed 
line in Figure 7. Within this region, assuming that void nucleation occurs, the size and shape of 
several cavitated voids at arbitrarily selected material points were determined by the deformation 
gradient The voids elongate in the direction of the applied loads and are close to the surface of the 
crack. 



^ applied ^ ^ 


Figure 7. A cavitation locus surrounding the deformed crack dp. 
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THE NONLINEAR ASYMPTOTIC SOLUTION FOR THE THIRD ORDER 

INVARIANT MATERIAL 


A second finite element analysis was completed using the third order invariant of the Rivlin 
constitutive model. Noting that Ij and I 2 are equal in plane strain, this material model can be 
expressed in terms of a single strain invariant, I, as 

W = (C 10 + C m ) (1-3) + (Cjo + C^ + Cjj) (/— 3) 2 + 

(£30 + ^03 + ^21 + ^ 12^ 

where Cy are the material constants. The nonlinear asymptotic solution 7 , under the assumptions 
of small scale nonlinearity, for this material model has been derived as 

— 324*/3r 7 ' / 6 2 7// 6 5—713 0 2 

y l 25 — (co(8) H- jCOse) x (-f(^;^;^;cos ej -^cos6) 


a f 1 — yr COS _ 1 /-2 

>2-^ /6 «in'[ l + „ m j (<o(e) + |cos0) - 

^^®r 7/ ®(co(8) + jcose) (1 + |sin 2 0-a>(8)cos0) 

, n - 0 (r- 2/3 ) 

t 22 -r -1 (?A) (|a) (<o(9) + ?COS0) 

T I2 - r _2/3 T 12 (0) + o(r -2/3 ) 

where 

G)(e) = ^ [9 -4sin 2 0] 1/2 
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THE NONLINEAR ASYMPTOTIC SOLUTION FOR THE THIRD ORDER 
INVARIANT MATERIAL (CONTINUED) 


2 2 1/2 
(1 + 2 sin 0 - <d( 6) cos0) 

cose — 

° 3^((o(0) + |cos0) 

A ~ C 21 + C 12 + C 30 + C 03 B = C ll + C 20 + C 02~ 9A 


and F[.„] and Tj2(0) arc hypergeometric functions. The sole load amplitude parameter, a, can 
again be determined by substituting the nonlinear asymptotic solution into the definition of the J- 
integral, yielding 


( 


a — 


3 <j) 

47t(lA 


2 \ 


1/6 




This asymptotic solution is remarkably different from that of the linear asymptotic solution and 
the nonlinear asymptotic solution for the Neo-Hookean material. The dominant stress, I s °f 
o(r !), but is inversely proportional to a function of 0. Both terms in the deformation field, y\ and 
y 2 , possess multiple functions of 0. 
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RESULTS FOR THE THIRD ORDER INVARIANT MATERIAL 


The finite element results will now be compared to the nonlinear asymptotic solution for this 
material. Circumferential plots for the y 2 - deformation component and the T 22 stress component 
surrounding the crack tip were examined in the refined mesh region, close to the crack tip. Again, 
the normalized results were given in reference to the undeformed coordinate system. Close agree- 
ment was found between the finite element results and the theoretical solution for this problem, 
as shown in Figures 8 and 9. 



Figure 8. Circumferential plot of y 2 taken along a radius of 4.6E-3 JappUed/F- 
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RESULTS FOR THE THIRD ORDER INVARIANT MATERIAL 

(CONTINUED) 


The circumferential plot of X22, as a function of 6, confirms the presence of the nonlinear 
asymptotic solution. A comparison of the finite element results to the linear elastic solution 
reveals some significant differences. In the linear elastic solution, the X22 stress component is 
higher in magnitude directly in front of the crack tip, when 6 is less than ninety degrees in the 
undeformed configuration. However for the nonlinear crack tip field, maxim u m values of T22 
were found as 6 approached 180 degrees in the undeformed configuration, close to the crack sur- 
face. Assuming that failure would occur at locations of maximum stress, these results suggest that 
material failure and secondary crack growth should initiate close to the crack surface, which is 
above and below the crack tip in the deformed configuration. The results were consistent with 
experimental observations on tearing in carbon black filled elastomers 8 . 



Figure 9. Circumferential plot of T 22 taken along a radius of 5.5E-3 Jappiied/M- 
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CONCLUSIONS 


A finite element analysis of a plane strain Mode I crack was performed for two different 
materials, a NeoHookean material and the third order invariant of the Rivlin constitutive law, 
under the assumptions of small scale nonlinearity. The numerical results for both ma t er ial s con- 
firmed the presence of the appropriate theoretical nonlinear asymptotic solution in a small region 
surrounding the crack tip. For the Neo-Hookean material, a potential cavitation locus encompass- 
ing the crack tip was identified. Within this cavitation region, possible sizes and orientation of 
nucleated microvoids were suggested. For the higher order material, localized material failure 
was predicted close to the crack surface, above and below the deformed crack tip, where the 
stress component reached maximum values. 
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INTRODUCTION 


Some recent developments in reduction techniques, as applied to predicting the tire contact response 
and evaluating the sensitivity coefficients of the different response quantities, are reviewed. The sensitivity 
coefficients measure the sensitivity of the contact response to variations in the geometric and material 
parameters of the tire. The tire is modeled using a two-dimensional laminated anisotropic shell theory with 
the effects of variation in geometric and material parameters, transverse shear deformation, and geometric 
nonlinearities included. The contact conditions are incorporated into the formulation by using a perturbed 
Lagrangian approach with the fundamental unknowns consisting of the stress resultants, the generalized 
displacements, and the Lagrange multipliers associated with the contact conditions. The elemental arrays are 
obtained by using a modified two-field, mixed variational principle. 

For the application of reduction techniques, the tire finite element model is partitioned into two regions. 
The first region consists of the nodes that are likely to come in contact with the pavement, and the second region 
includes all the remaining nodes. The reduction technique is used to significantly reduce the degrees of 
freedom in the second region. 

The effectiveness of the computational procedure is demonstrated by a numerical example of the 
frictionless contact response of the space shuttle nose-gear tire, inflated and pressed against a rigid flat surface. 
Also, the research topics which have high potential for enhancing the effectiveness of reduction techniques 
are outlined. 
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OBJECTIVES AND SCOPE 


Despite the significant advances made in the development of effective computational models and 
computational strategies for the numerical simulation of the nonlinear tire response (see, for example, Refs. 
1 , 2 and 3), the modeling of aircraft tire contact response with the pavement remains one of the challenging 
applications of computational structural mechanics. The objectives of this presentation are to (see Fig. 1): a) 
review some of the recent developments in reduction techniques as applied to the tire contact problem, and 
b) identify some of the future directions for research. The scope of the presentation covers both the frictionless 
contact response of tires as well as sensitivity analysis. 


Objectives 

• Review some of the recent developments in reduction techniques 
as applied to the tire contact problem 

• Identify some of the future directions for research 
Scope 

• Frictionless contact response of tires 

• Sensitivity analysis 
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FINITE ELEMENT MODELING 
Two-Dimensional First-Order Shear Deformation Theory 


The analytical formulation for frictionless contact of tires is based on a form of the nonlinear, moderate 
rotation, Sanders-Budiansky type shell theory (Refs. 4 and 5) with the effects of variations in the geometric 
and material parameters, anisotropic material response, and transverse shear deformation included. The 
fundamental unknowns consist of the stress resultants, the generalized displacements, and the Lagrange 
multipliers associated with the contact conditions (see Fig. 2). 


h 



Figure 2. 
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FINITE ELEMENT MODELING (Cont'd.) 


The discrete finite element equations governing the frictionless contact response of the tire are obtained 
by using a modified, two-field, Hellinger-Reissner mixed variational principle. The modification consists of 
augmenting the functional of that principle by two terms: the Lagrange multiplier associated with the nodal 
contact pressures and a regularization term which is quadratic in the Lagrange multipliers (Refs. 6 to 9). The 
resulting nonlinear equations are shown in Fig. 3. The vector <Z> is the global response vector which includes 
the stress-resultant parameters and the nodal values of the generalized displacements; {A > is the vector of nodal 
values of Lagrange multipliers associated with the contact nodes; [K] is the linear matrix of the discrete tire 
model; [Q], [R] are the matrices associated with the contact condition and the regularization term; <G(Z)} is 
the vector of nonlinear terms; p is a load parameter; <P> is the vector of normalized external forces; <g 0 > is 
the vector associated with initial gaps in the contact region; and £ is a penalty parameter. 


Response (Nonlinear Equations) 
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Figure 3. 
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FINITE ELEMENT MODELING (Cont'd.) 


The governing equations for the sensitivity coefficients are shown in Fig. 4. In these equations, d 
represents a typical tire parameter. Note that the governing equations for the response are nonlinear but those 
for the sensitivity coefficients are linear. 


Sensitivity Coefficients (Linear Equations) 




Figure 4. 


SENSITIVITY COEFFICIENTS 


Figure 5 highlights the sensitivity coefficients for tires, their uses, and the approaches used for their 
evaluation. 

The tire parameters, with respect to which the sensitivity of the response need to be evaluated, include 
rubber and cord material properties, cord end counts, epi, cord diameters, and cord angles with the meridional 
direction. 

Among the uses of the sensitivity coefficients are: 

1 . Determine a search direction in the application of nonlinear mathematical programming algorithms. 

2. Assess effects of uncertainties in geometric and material parameters of the tire. 

3. Predict changes in tire response due to changes in material parameters. 

4. Identify changes in material and geometric parameters required to achieve certain performance 
characteristics. 

The different approaches used in evaluating the sensitivity coefficients can be grouped into the following 
three categories: 

1 . Analytical direct differentiation methods. 

2. Finite different methods. 

3. Semi-analytical or quasi-analytical methods. 

Tire Parameters 

d refers to rubber and cord material properties, epi, cord diameters and 
cord angles 

Uses 

• Determine a search direction in the application of nonlinear mathematical 
programming algorithms 

• Assess effects of uncertainties in geometric and material parameters of 
the tire 

• Predict changes in tire response due to 
changes in material parameters 

• Identify changes in material and 
geometric parameters required to 
achieve certain performance 
characteristics 

A pproaches 

• Analytical direct differentiation methods 

• Finite difference methods 

• Semi-analytical or quasi-analytical methods 



Figure 5. 
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REDUCTION TECHNIQUE FOR EVALUATING THE RESPONSE 
AND SENSITIVITY COEFFICIENTS 


The reduction technique for evaluating the response and sensitivity coefficients is highlighted in Fig. 6. 

The response vector {Z} is approximated by a linear combination of a few preselected basis vectors 
(columns of the matrix [T]). The approximation can be thought of as transformation of variables (from {Z} 
to {t|/}). The first-order derivatives of {Z} with respect to d are obtained by differentiating the approximate 
expression of {Z} with respect to d. 

An effective choice for the basis vectors was found to be the path derivatives; namely, the various-order 
derivatives of {Z> with respect to a control parameter p. 


Basis Reduction 


az\ -fri 
ad/ ni1 - [ MadJ rj1 


r « n 


ar 


ad 


n,r 


Wr,1 


where [r] = transformation matrix 

{i| /}, = unknown parameters 


Selection of m 



p = load (or displacement parameter) 


S = o 



Figure 6. 
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REDUCTION TECHNIQUE FOR EVALUATING THE RESPONSE 
AND SENSITIVITY COEFFICIENTS (Cont’d.) 

Treatment of Contact Response 


For the application of the reduction technique to the tire contact problem, the finite element model is 
partitioned into two regions. The first region consists of the contact region (i.e., the nodes that are likely to 
come in contact with the pavement) and the second region consists of all the remaining nodes (see Fig. 7). The 
response vector {Z} is partitioned accordingly. 

The subvector {Z { } is expected to exhibit a significant change with the change in the applied load, and 
therefore, is not approximated. By contrast, the subvector {Z 2 } is not expected to change much, and is 
approximated by a linear approximation of a few preselected vectors (columns of the matrix [T]). 

A Bubnov-Galerkin technique is used to replace the governing equations for the contact response and 
its sensitivity coefficients by two reduced systems of equations: a nonlinear system of equations in {Zj }, {\jr}; 

/3Z, \ /chjr\ 

and a linear systems of equations in V -gj-j, • 
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Figure 7. 
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SPACE-SHUTTLE NOSE-GEAR TIRE 
Characteristics of the Cross Section 


To test and evaluate the effectiveness of the reduction technique, the frictionless contact response of the 
space shuttle orbiter nose-gear tire and its sensitivity coefficients were generated by this technique. The 
geometric characteristics of the tire cross section are shown in Fig. 8. Each half of the tire cross section was 
divided into seven segments. Each segment contained a different number of layers, different material 
properties corresponding to different cord content in the composite, and varying cord orientations. 



Figure 8. 
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STATIC LOAD DEFLECTION CHARACTERISTICS 


Figure 9 shows a plot of the applied displacement versus the total vertical contact load obtained by two 
different discrete models. The full inflation pressure of 320 psi was applied to the tire, then the inflated tire 
was pressed against the rigid pavement, and the applied displacement was increased. Model 1 has a total of 
10,188 displacement degrees of freedom in the full tire and model 2 has a total of 20,736 displacement degrees 
of freedom. The close agreement between the predictions of the two finite element models and the experiments 
is clearly seen in Fig. 9. 
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Figure 9. 



ACCURACY OF NONLINEAR TIRE RESPONSE PREDICTED BY REDUCTION TECHNIQUE 


Figure 1 0 shows the accuracy of the total vertical contact load and the total strain energy U, in the region 
-0.45 < ^ < 0.45 obtained by using the reduction technique in conjunction with model 2. The tire was divided 
into two regions. The contact region (region 1) had 432 elements and 9065 displacement degrees of freedom. 
Region 2 had 648 elements and 1 1,31 1 displacement degrees of freedom. The basis vectors were generated 
at A/h 0 = 0.1 and were not updated throughout the range of the applied displacement considered. For the range 
of A/ho considered, the total vertical contact force obtained by using ten basis vectors are almost undistin- 
guished from those obtained by using the full system of finite element equations. The error in the total strain 
energy obtained by using ten vectors was less than 2% (the scale selected in Fig. 10 amplifies that error). 



Figure 10. 
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VARIATION OF TOTAL STRAIN ENERGY DENSITY DISTRIBUTION 
WITH APPLIED DISPLACEMENT 


Figure 1 1 shows the variation of the strain energy density with applied displacement in the region -0.45 
< ^ < 0.45. The strain energy density contours are normalized by dividing them by Ex ho- The strain energy 
density is primarily influenced by the inflation pressure. Since the normal tractions associated with contact 
are compressive, the strain energy density is reduced in the contact region as well as in the sidewall area near 
the bead (away from the contact region). 



Figure 1 1 . 
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ACCURACY OF THE SENSITIVITY COEFFICIENTS OF THE TOTAL STRAIN ENERGY 
WITH RESPECT TO CORD DIAMETERS AND YOUNG'S MODULI 


Figure 1 2 gives an indication of the accuracy of the total strain energy with respect to the cord diameters 
d] and d 2 , and Young's moduli of the cord and rubber E c and E r . As can be seen from Fig. 12, reasonable 
accuracy is obtained with only ten basis vectors. 



Figure 12. 




ACCURACY OF THE SENSITIVITY COEFFICIENTS OF THE TOTAL STRAIN ENERGY 
WITH RESPECT TO THE PARAMETER OF THE CORD END COUNT AND CORD ANGLE 


Figure 1 3 gives an indication of the accuracy of the sensitivity coefficients of the total strain energy with 
respect to the parameters of the cord end count, epi, and the cord angle 0 . Again, for the entire range of applied 
displacements, reasonable accuracy is obtained with only ten vectors. 
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Figure 13. 
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NORMALIZED CONTOUR PLOTS OF TOTAL STRAIN ENERGY DENSITY 


au /au 

8d / 5d 


max 


Figure 1 4 shows normalized contour plots of the sensitivity coefficients of the total strain energy density 
U with respect to four tire parameters (E c , E r , d j and d 2 ) at A/h 0 =0.8 and A/ho =3.2. Each sensitivity coefficient 
is normalized by dividing by the maximum value at the given A/h 0 . The distributions of the sensitivity 
coefficients of U differ from each other and from the corresponding U distributions. Also, the distributions 
of the sensitivity coefficients away from the contact area remain nearly axisymmetric. 
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Figure 14. 













FUTURE DIRECTIONS OF RESEARCH ON REDUCTION TECHNIQUES 


Some of the research topics which have high potential for enhancing the effectiveness of reduction 
techniques are (see Fig. 15): 

1 . Development of reduction techniques for hierarchical sensitivity analysis, which includes evaluation 
of the sensitivity coefficients of the tire response with respect to a hierarchy of parameters ranging 
from the stiffnesses of the two-dimensional laminated shell model of the tire to the material and 
geometric parameters of the cord, rubber and interface/interphase. 

Hierarchical sensitivity analysis can be used to: 

a) bridge the gap between material development and structural design of tires; 

b) study the response of tires at different length scales; and 

c) study the effect of imperfections in the tire materials, cord diameters, and cord angles on the 
response. 

2. Application of reduction techniques to: a) damage and failure analysis; b) thermomechanical 
response of tires; and c) rolling tire problem. 


• Development of reduction techniques for hierarchical 
sensitivity analysis 

- Bridge gap between material development and structural 
design of tires 

- Study the response of tires at different length scales 

- Study of the effect of imperfections on the tire response 

• Applications to 

- Damage and failure analysis 

- Thermomechanical response of tires 

- Rolling tire problem 

Figure 15. 
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SUMMARY 


In summary, some recent advances in the application of reduction techniques to tire contact problems 
have been reviewed. This included both prediction of contact response and sensitivity analysis. 

Reduction techniques have the potential of making detailed (large scale) contact and sensitivity analysis 
practical. 

Among the research topics that have high potential of increasing the effectiveness of reduction tech- 
niques are: a) hierarchical sensitivity analysis, which can be used to bridge the gap between material 
development and structural design of tires and to help in identifying the length scale required in modeling the 
different regions of the tire; and b) application of reduction techniques to damage and failure analysis of tires. 


• Reviewed some recent advances in the application of 
reduction techniques to tire contact problems 

- prediction of contact response 

- sensitivity analysis 

• Reduction techniques can make detailed (large scale) 
contact and sensitivity analyses practical 

• Future directions of research on reduction techniques 

- hierarchical sensitivity analysis-bridge gap between 
material development and structural design of tires 

- application to damage and failure analysis of tires 


Figure 16. 
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TIRE FOOTPRINT STUDIES 


Mangal Chawla and John Medzorian 
Landing Gear Systems Group 
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Wright Patterson Air Force Base, OH 


ABSTRACT 


This presentation covers the results of tire footprint studies conducted in the Landing Gear Develop- 
ment Facility of the USAF Wright Laboratory at the Wright Patterson Air Force Base, OH. Tire footprint 
studies are essential in understanding tire wear mechanisms and computing tire tread wear rates. The power 
input into the tread is the driving force for tread wear. Variables needed for power input calculations include 
the footprint pressure and slip velocity distributions. Studies were performed on the effects of power input 
distributions due to vertical load, camber, yaw, inflation pressure and tire construction. For the present study, 
two tire constructions, one radial and the other bias, were selected. These tires were for the F-16 Block 30 
Fighter aircraft, both of which were previously worn. The present study was limited to steady straight roll with 
a 14,000 lb vertical load, a 310 psi inflation pressure, zero yaw and camber. All tests were conducted on the 
Tire Force Machine (TFM) with a specialized sensor plate with embedded pressure sensors (X, Y and Z) and 
slip sensors (X & Y). All tests were conducted for atable speed of 1 in/s. Tests on the TFM show that the power 
intensity distributions and total power for both tire constructions are quite similar for straight roll. Later on, 
tests were also conducted on a modified dynamometer which was overlaid with a grit wear surface. The tire 
speed was maintained at 40 miles per hour and yaw was set to four degrees. Dynamometer tests showed that 
radial tires have more tread wear than the bias tire, however, in the field, radial tires have longer life. This is 
a dilemma. 
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PROBLEM 


Included is the problem statement of low tire life. 

■ Excessive tire wear due to 

- Large contact pressures 

- Large slip velocity 

■ Very few sorties/tire 

- F-16 Block 30.. . 8- 12 landings 

- F-16 Block 40 & 50.. . 4-8 landings 

Figure 1 . 


92 



OBJECTIVES 


This slide states the objectives of the program. The main objectives are to quantify parameters affecting 
tread wear. 


■ Overall Objective: Increase tread life 

u Focussed Objective: Quantify operating conditions and tire design 
parameters affecting tread wear 

• Operating conditions 

- Load 

- Takeoff/Landing profile (speed, yaw, camber, braking) 

- Antiskid characteristics 

• Tire design parameters 

- Inflation pressure 

- Tire construction (bias/radial) 


Figure 2. 



NEEDS 


This slide states the need to study footprint pressures and slips so that power intensity can be calculated 
tread-wear computations. This data is needed for steady roll, yaw and braking conditions. 


■ Determine aircraft tire footprint slips and pressures 
to quantify power intensity due to 

- Steady roll 

- Yaw 

• Steady (crosswind) 

• Cornering 

- Braking 

• Steady (constant slip) 

• Transient (antiskid) 


Figure 3. 



TOOLS 


The tools for carrying out this work are the Tire Force Machine (TFM) and specialized sensors which 
measure footprint pressures and slips. 

■ Experimental 

- TFM 

• Steady conditions only 

• Low speed (2 in/s max) 

m Analytical/Computational 

- FEM 

Figure 4. 
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APPROACH 


The approach to this work is to perform roll tests on the TFM and measure footprint data using specialized 
equipment. The data will be analyzed for wear energy rate determination. Parameter studies will be performed 
to compare the performance of radial against bias tire for vertical load, yaw and camber. 


■ Perform Footprint measurements on TFM 

- Slip velocity x (longitudinal) y (lateral) 

- Contact force intensity x, y and z (normal) 
m Data Analysis 

- Wear energy rate determination 

- Parametric studies on F-16 tires 

* Type (Radial vs Bias) 

* Vertical Load 

* Yaw 

* Camber 


Figure 5. 



SENSORS FOR TIRE FORCE MACHINE 


This slide shows the layout of pressure and slip sensors and the direction of roll. Pressure sensors give 
X, Y and Z components of pressure, whereas slip sensors measure X and Y components. There are 37 pressure 
gages and 37 slip gages. The holes are 1/4 inch in diameter. The pressure sensors are level with the top of the 
plate whereas the slip sensors are 4 mils above the plate surface. 


SENSORS FOR TIRE FORCE MACHINE 

I .SO' X 36 SPACES - 18' 



Figure 6. 
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F-16 BLOCK 30 


This is a mesh plot of the vertical (Z) pressure for the bias and the radial F-16 Block 30 tires. Both plots 
look similar in shape and magnitude indicating that the normal pressures are nearly the same for these tires. 
In all runs, the tire is rolling towards the origin in the fore-aft direction. The table speed was 0.98 1 in/s for the 
radial tire and 0.913 in/s for the bias tire. 


F-16 BLOCK 30 

MAS Load = 14 Kips, Inf. Pr. = 310 PSI RADIAL 


Yaw = 0, Camber = 0 


25.5x8.0-14/20 load=14.0001b P=310psi Y=C=0 


25.5x8.0R14/20 Load=14,000lb P=310psi Y=C=0 



(ore-aft position (in) 


0 0 


lateral position (in) 


fore-aft position (in) 


lateral position (in) 


Figure 7. 


X PRESSURE 


This is a plot of the fore-aft (X) pressure distributions for both radial and bias tires. The total measured 
steady state X force is 1 40 lbs for the radial tire and 271 lbs for the bias tires. This compares favorably with 
the predicted integrated values of 265 lbs and 1 41 lbs. Therefore, the total dissipative power in these runs works 
out to be 137 in-lb/s for the radial tire and 247 in-lb/s for the bias tire. 


F-16 BLOCK 30 

BIAS j_ oad _ 14 Kips, )nf p r _ 310 pgj RADIAL 


Yaw = 0, Camber = 0 



Figure 8. 
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Y pressure (psi) 


Y PRESSURE 


This is a plot of the lateral (Y) pressure distribution for both radial and bias tires. 

F-16 BLOCK 30 

BIAS Load = 14 Kips, Inf. Pr. =310 psi RADIAL 
Yaw = 0, Camber = 0 


25.5x8.0-14/20 Load=14,000lb P=310psi Y=C=0 25.5x8.0R14/20 Load=14,000lb P=310psi Y=C=0 



fore-aft position (in) 0 0 )atera i posjt j on (j n ) fore-aft position (in) 0 0 lateral position (in) 


Figure 9. 
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TRACTION VECTOR 


This is a plot of the traction vectors providing both magnitude and direction. 


F’16 BLOCK 30 

BIAS Load = 14 Kips, Inf. Pr. = 310 PSI RADIAL 
Yaw = 0, Camber = 0 


25.5x8.0-14/20 Load=14,000lb P=3t0psl Y=C=0Traction vector 



25.5x8.OR14/20 Load-14, OOfflb P=310psi Y=C=0Traction vector 



Figure 10 





X SLIP VELOCITY 


This is the slip velocity in the fore-aft (X) direction, measured in mil/s. Note that the peaks occur at the 
trailing edges of the contact patch. The slip velocity was calculated from the slips by using the central 
difference method. The data was manipulated such that at zero normal (Z) pressure, no slip velocity occurs. 
Note that the maximum slip occurs at the trailing edge of the tire. 

F-16 BLOCK 30 

BIAS Load = 14 Kips, Inf. Pr. = 310 PSI RADIAL 


Yaw = 0, Camber = 0 

25 5x8.0-14/20 Load=14.000lb P=310p*i Y=C=0 25.5x8 0R14/2O LoM=14,OOOIb P-SIOpsi Y=C=0 



(ore-aft position (in) 0 0 latera) position (in) fore-aft position (in) 0 0 Jatefa( posjti()n ^ jn j 


Figure 1 1 . 
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Y slip velocity mil/s 


Y SLIP VELOCITY 


This is the slip velocity in the lateral (Y) direction, measured in mil/s. Note that the maximum values 
occur at the trailing and leading edges of the contact patch. 


F-16 BLOCK 30 

BIAS Load = 14 Kips, Inf. Pr. =310 psi RADIAL 
Yaw = 0, Camber = 0 



fore-aft position (in) 0 0 lateral position (in) fore-aft position (in) 0 0 lateral position (in) 


Figure 12. 
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slip velocity magnitude (mil/s) 


SLIP VELOCITY MAGNITUDE 


This is the slip velocity magnitude measured in mil/s. Note that the peaks occur at the trailing edge of 
the contact patch. 


F-16 BLOCK 30 

BIAS |_ oad _ .J4 Kj p S| | nf Pr _ 3 -)o psi RADIAL 


Yaw = 0, Camber = 0 

25.5x8.0-14/20 Load=14,000lb P=310psi Y=C=0 25.5x8.0R14/20 Load=14,000lb P=310psi Y=C=0 



Figure 13. 


SLIP VECTOR 


This is the plot of the slip velocity vector giving both the magnitude and the direction. There is a lot of 
similarity between the two tires. 


F-16 BLOCK 30 

BIAS Load = 14 Kips, Inf. Pr. =310 PSI RADIAL 
Yaw = 0, Camber = 0 

25.5x8.0-1420 l_oad=14, 0001b Yaw=Cam=0 Slip vector 25,5x8.0R14/20 Load=14,000lb P=310psl Y=C=0 Slip vector 



Figure 14. 
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Z PRESSURE 


This is a contour plot of the normal (Z) pressure for the bias tire. When integrated over the contact area, 
it gives a total value of 15,565 lbs which compares favorably with the applied load of 14,000 lbs. 

F-16 BLOCK 30 BIAS 

Load = 14 Kips, Inf. Pr. = 310 PSI 
Yaw = 0, Camber = 0 

Pressure 
(psi) 

1 483.556 
423.111 
362.667 
302.222 
241.778 
181.333 
W 120.889 
60.4444 


2 4 6 8 10 12 

Fore-aft position (in) 


Figure 15. 
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SLIP VELOCITY MAGNITUDE 


This is a contour plot of the slip velocity magnitude for the bias tire. Note that the slip velocity is 
maximum at the trailing edge of the contact patch. This data was manipulated by zeroing out values which 
occurred at zero contact pressure. 


Slip Vel Mag 
(mil/s) 

I 62.3822 
I 46.3498 
I 34.4378 
I 25.5872 
I 19.0112 
Wj 14.1253 
W 10.495 
7.79778 


2 4 6 8 10 12 

Fore-aft position (in) 


Figure 16. 


F-16 BLOCK 30 BIAS 


Load = 14 Kips, Inf. Pr. =310 PSI 
Yaw = 0, Camber = 0 
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POWER INTENSITY 


This is a contour plot of the power intensity calculated by multiplying an assumed value of coefficient 
of friction (0.5) and the Z pressure by the slip velocity magnitude. Note that the maximum power intensity 
occurs at the front and back edges of the tire. The total dissipative power due to friction was determined by 
integrating this power intensity over the contact area. The resulting value was 22.0 in-lb/s for the bias tire. 


Pow Int 
i-lb)/(s in A 2) 

3.71644 
3.25189 
2.78733 
2.32278 
1 .85822 
1 .39367 
0.9291 1 1 
0.464556 


2 4 6 8 10 12 

Fore-aft position (in) 


Figure 17. 


F-16 BLOCK 30 BIAS 

Load = 14 Kips, Inf. Pr. =310 PSI 
Yaw = 0, Camber = 0 
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Z PRESSURE 


This is a contour plot of the normal (Z) pressure for the radial tire. When integrated over the contact area 
it gives a value of 15,095 lbs which compares favorably with the applied load of 14,000 lbs. 


F-16 BLOCK 30 RADIAL 

Load = 14 Kips, Inf. Pr. = 310 PSI 
Yaw = 0, Camber = 0 

Pressure 
(psi) 

■ 464.889 
I 406.778 
I 348.667 
I 290.556 
H 232.444 
|| 174.333 
116 222 
|“j 58.1111 


Fore-aft oosition 8 fin) 10 


Figure 18. 
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SLIP VELOCITY MAGNITUDE 


This is a contour plot of the slip velocity magnitude for the radial tire. Note that the slip velocity is 
maximum on the trailing edge of the contact patch. This data was manipulated by zeroing out values which 
occurred at zero contact pressure. 


F-16 BLOCK 30 RADIAL 


Load = 14 Kips, Inf. Pr. =310 PSI 
Yaw = 0, Camber = 0 


Slip Vel Mag 



Fore-aft position 8 (in) 


Figure 19. 
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POWER INTENSITY 


This is a contour plot of the power intensity for the radial tire. This was calculated by multiplying an 
assumed value of coefficient of friction (0.50) by the Z pressure and the slip velocity magnitude. Note that 
the maximum power intensity occurs at the front and back edges of the tire. The total dissipative power was 
23.6 in-lb/s for the radial tire. 


F-16 BLOCK 30 RADIAL 

Load = 14 Kips, Inf. Pr. =310 PSI 


Yaw = 0, Camber = 0 



Pow Int 
(in-lb)/(s in A 2) 


4.13244 

3.07039 

2.28129 

1.69499 

1.25938 

0.935713 

0.695233 

0.516556 


Fore-aft position 8 (in) 


Figure 20. 


ill 




POWER INTENSITY 


This is a mesh plot of the power intensity distributions for the two tires. 


F-16 BLOCK 30 

= 14 Kips, Inf. Pr. = 3 
Yaw = 0, Camber = 0 


BIAS Load = 14 Kips, Inf. Pr. =310 PSI RADIAL 


25.5x8.0-14/20 Load=14,000lb P=310psi Y=C=0 



25.5x8.0R1 4/20 Load=14,000lb P=310psi Y=C=0 


fore-aft position (in) 


lateral position (in) 



Figure 21. 
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POWER PER UNIT LENGTH 

This is a plot of power per unit width plotted in the lateral direction which was determined by integrating 
the power intensity along the fore-aft direction. 

F-16 BLOCK 30 

BIAS Load = 14 Kips, inf. Pr. = 310 PSI RADIAL 
Yaw = 0, Camber = 0 


25.5*6.0-14/20 Load*14,000«> P*310pii Y«C*0 25.5xa.0R14/20 Load=14.000U> P«310psi Y<C>0 



UUcXpoUion(n) laMral position (In) 


Figure 22. 
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RESULTS 


Tabulated results for both bias and radial tires. All values are very similar. 



BIAS 

RADIAL 

Table speed (in/s) 

.913 

.981 

Maximum slip (mil) 

31.3 

49.1 

X-Pressure range psi 

(-122,94) 

(-147,106) 

Y-Pressure range psi 

(-113,126) 

(-132,115) 

Traction Vector psi 

(0,129) 

(0,148) 

Z- Pres sure range psi 

(0,544) 

(0,523) 

X-Slip Velocity mil/s 

(-70,50) 

(-63,64) 

Y-Slip Velocity mil/s 

(-25,30) 

(-35,35) 

Slip Velocity Magnitude mil/s 

(0,70.2) 

(0,66.1) 

Power Intensity in-lb/(s-in A 2) 

(0,4.2) 

(0,4.6) 


MEASURED (INTEGRATED) TOTAL VALUES 


Fx (lb) 

271 (265) 

140 (141) 

Fy (lb) 

182 (6) 

157 (199) 

Fz (lb) 

14,000 (15,565) 

14,000 (15,291) 

Power (in-lb/s) 

(22.0) 

(23.6) 


Figure 23. 



WL DYNAMOMETER TIRE WEAR RESULTS SUMMARY 


This is a plot of the actual tire wear on the bias and the radial tires which were run on a modified external 
drum dynamometer overlaid with an abrasive material glued to the wheel surface. The tire was run with four 
degrees of yaw at 40 mph. Tread wear was quantified by measuring the maximum groove depth on the tire 
tread. Note that on the dynamometer for four degrees of yaw, the radial tire wears much faster than the bias 
tire. However, in the field, the radial tires have longer tread life. 



Figure 24. 





Determination of Tire Cross-Sectional Geometric 
Characteristics From a Digitally Scanned Image 

Kent T. Danielson 

Center for Computational Structures Technology 
University of Virginia 
Hampton, VA 


117 




DETERMINATION OF TIRE CROSS-SECTIONAL GEOMETRIC 
CHARACTERISTICS FROM A DIGITALLY SCANNED IMAGE 


Kent T. Danielson 

Center for Computational Structures Technology 
University of Virginia 
NASA Langley Research Center 
Hampton, Virginia 23681 


INTRODUCTION 


A semi-automated procedure is described for the accurate determination of geometrical 
characteristics using a scanned image of the tire cross-section. The procedure can be useful for 
cases when CAD drawings are not available or when a description of the actual cured tire is 
desired. 

The determination of geometric details from digitally scanned images has been used for various 
applications (e.g., air flow behavior [1]). It is believed, however, that this is the first published 
application to tires. Curves representing the perimeter of the tire cross-section are determined by 
an edge tracing scheme, and the plyline and cord-end positions are determined by locations of color 
intensities. Edge and contour tracing techniques and their resultant chain codes were introduced by 
Freeman [2]. Since then, numerous variations have been reported. Several survey articles have 
reviewed and analyzed these techniques [3-5], 

The current procedure provides an accurate description of the perimeter of the tire cross-section 
and the locations of plylines and cord-ends. The position, normals, and curvatures of the cross- 
sectional surface are included in this description. The locations of the plylines provide the 
necessary information for determining the ply thicknesses and relative position to a reference 
surface. Finally, the locations of the cord-ends provide a means to calculate the cord-ends per inch 
(epi). 

Menu driven software has been developed to facilitate the procedure using the commercial code, 
PV-Wave by Visual Numerics, Inc., to display the images. From a single user interface, separate 
modules are executed for image enhancement, curve fitting the edge trace of the cross-sectional 
perimeter, and determining the plyline and cord-end locations. The code can run on SUN or SGI 
workstations and requires the use of a mouse to specify options or identify items on the scanned 
image. 
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OBJECTIVES 


The objective is to produce a robust procedure with the following characteristics: 


1 . An accurate determination of the 
cross-sectional geometry of an actual tire 

a) Perimeter shape 

b) Ply centerline locations 

c) Ply end counts 

2. Simple user-friendly software 

3 . Produce a simple database that can be 
used to create finite element models, to 
merge with CAD systems, or to be 
manipulated for other purposes. 
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RAW SCANNED IMAGE 


The procedure begins by capturing the tire cross-sectional image in monochrome with a flat bed 
scanner. A raw scanned image of the space shuttle nose gear tire is shown in Figure 1 . The image 
is a grid of points (pixels), each with its own grayscale value (scaled numerical value specifying 
amount of black and white). The raw image is digitally enhanced for use with the software. A 
slight enhancement is used to more easily inspect the accuracy of the calculated geometry with 
respect to the tire image. In addition, a higher degree of enhancement is used with the algorithms 
that calculate the geometric quantities. 



Figure 1 
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DIGITALLY ENHANCED IMAGES 


A separate module is executed from the user interface to produce enhanced images. The raw 
image is digitally enhanced by subtracting a smooth area-averaged copy, which has a 
predominance of low frequency components of the image, from the original. By subtracting out a 
smoothed or averaged image, the high frequency components, especially the tire cord-end images, 
are enhanced. A highly enhanced image corresponding to Figure 1 is shown in Figure 2. The 
highly enhanced image is used for the calculation of the geometric quantities. A slightly enhanced 
image will also be used to assess accuracy by superimposing the calculated geometry on it. The 
degree of enhancement is determined by the size of the averaged area. The software allows the 
user to specify the parameters that determine the degree of enhancement, and the enhanced images 
are displayed for inspection. 



Figure 2 
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DETERMINATION OF CROSS-SECTIONAL PERIMETER 


A separate module is also executed from the user interface to fit splines to the perimeter of the 
tire cross-section. The grayscale values are used to determine geometrical details, since geometry 
changes will generally coincide with changes in color. An edge tracing procedure is used to 
determine the inner and outer surface by determining when pixel grayscale values go from white to 
black. After the edge trace is completed, B-splines are fit to the inside and outside surfaces. The 
curves are used to calculate coordinates, normals to the surface, and curvatures. The software 
allows the user to change the spline order, number of splines, and spacing of splines. For 
inspection, the splines are then graphically superimposed over a plot of the edge trace points. 
Figure 3 shows the B-splines fit to the edge trace of the space shuttle nose gear tire. 
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DETERMINATION OF PLY CENTERLINES 


A separate module is also executed from the user interface to determine the locations of the ply 
centerlines and cord-ends. Similar to the edge trace procedure, the tire cord locations within the 
rubber will be determined by grayscale summations that indicate a cluster of white in the image. 
The plyline locations for different positions along the tire meridian are determined by a series of 
points on the perimeter. A mouse click on each desired point on the outside surface is made on the 
image, and the normal to the outside surface is determined at the point. Figure 4 shows an 
example for point A, which is also identified in Figure 3. The distance to the inside surface along 
the normal is calculated, and the distance to each ply centerline is determined by a summation of 
pixel grayscale values. Since the normal will not generally cross all cord-ends, a sum of pixel 
values is made for a fixed distance on each side of the normal. For the most accurate sums, the 
summation direction on each side of the normal would ideally be along the direction of the ply. 
The code provides a simple, yet effective, way of specifying this direction. The option is given for 
using the outside surface, the inside surface, or a smooth transition between the two. 



Figure 4 
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Grayscale Value Sum, 1 x e3 


DETERMINATION OF PLY CENTERLINES 


Each grayscale sum is assigned to its pixel position along the normal. When plotted along the 
normal, as shown in Figure 5 for point A, these sums should show peaks of white where the 
plylines lie. The peaks are chosen to be the ply centerlines and then plotted on the image for 
inspection. The first two peaks in Figure 5 have smaller amplitudes than the other peaks and are a 
result of noise in the heavily enhanced image. Inaccurate plylines like these can be deleted or new 
desired ones added with simple mouse clicks on the image. The other ten peaks represent the ply 
centerlines. Figure 4 shows the ply centerlines, as determined from the ten major peaks in Figure 
5, marked by crosses on a slightly enhanced image. This procedure is repeated for other points 
along the tire meridian, until an adequate description is obtained. 



Figure 5 



DETERMINATION OF CORD-END SPACINGS 


In a similar manner, cord-end spacings can also be obtained for the plies determined at each of 
the points along the tire meridian. Since the plies are oriented at an angle to the cross-section, 
projections of the actual spacings are calculated. At each of the ply centerlines, the pixel grayscale 
values are plotted for the fixed distance on both sides of the normal. The peaks are chosen as the 
cord-ends and plotted on the image. Inaccurate cord-ends can again be deleted or new ones added 
with mouse clicks on the image. Since the cord-ends will be calculated to lie on a given line, an 
option to choose the line is provided. Besides the inner surface, the outer surface, or a smooth 
transition between the two, a straight line may alternatively be specified as the ply direction by 
making mouse clicks on each side of the normal shown on the image. The cord-end spacing is 
determined as an average distance between cord-ends. Figure 6 shows four cord-ends for the 
second ply of Figure 5, as determined from peaks in the grayscale values, marked by crosses on a 
slightly enhanced image. 



Figure 6 
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SUMMARY 


A semi-automated procedure was described for the accurate determination of geometrical 
characteristics using a scanned image of the tire cross-section. The procedure is useful for cases 
when CAD drawings are not available or when a description of the actual cured tire is desired. In 
general, the automated features are not entirely adequate and the program will therefore require 
some user intervention. Although accurate cross-sectional details are automatically determined 
most of the time, the ultimate control of each detail is given to the user. For efficient use of the 
automation algorithms, the user also has continuous access to the parameters that control them. 
The calculated geometry is plotted on the image to ensure accuracy. If necessary, details can then 
be modified with simple options available in the software. The software produces the following: 

1 . An accurate determination of the 
cross-sectional geometry of an actual tire 

a) B-splines for the perimeter shape 

b) Ply centerline locations 

c) Ply end counts 

2. Simple user-friendly software using 
PV-Wave to display the images 

a) Menu driven 

b) Single user interface to execute 
separate modules 

c) Graphical display of calculated details 
on scanned image 

3. A simple database that can be used to 
create finite element models, to merge 
with CAD systems, or to be manipulated 
for other purposes. 
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INTRODUCTION 


The modeling of tires and the simulation of tire behavior are complex problems. The 
MacNeal-Schwendler Corporation (MSC) has a number of finite element analysis 
products that can be used to address the complexities of tire modeling and simulation. 
While there are many similarities between the products, each product has a number of 
capabilities that uniquely enable it to be used for a specific aspect of tire behavior. 

• MSC/NASTRAN - General purpose finite element program for 
linear and nonlinear static and dynamic analysis. 

• MSC/ABAQUS - Nonlinear statics and dynamics finite element 
program. 

• MSC/PATRAN AFEA (Advanced Finite Element Analysis) - 
General purpose finite element program with a subset of linear and 
nonlinear static and dynamic analysis capabilities with an 
integrated version of MSC/PATRAN for pre- and post-processing. 

• MSC/DYTRAN - Nonlinear explicit transient dynamics finite 
element program. 
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TIRE CHARACTERISTICS AND GEOMETRY 


Tires are very complex structures that have fairly simple overall geometry, but have 
complicated cross-sections (Figure 1). In their operating environment tires are subjected 
to varying types and magnitudes of loads. Some of the important aspects of tire behavior 
can be summarized as: 

• Tire Footprint and Contact Area 

• Dynamics (landing, curb impact, etc.) 

• Cornering Loads 

• Vibration 

• Hydroplaning 

• Interaction with Road Surface 




Figure 1. Typical Tire Geometry and Cross Section Detail 
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IMPORTANT ASPECTS IN SIMULATING TIRE BEHAVIOR 


There are many important aspects of tire behavior which must be considered in a tire 
simualtion: 


• Contact and Friction 

• Material Constitutive Models (Linear and Nonlinear Behavior) 

• Geometric Nonlinear Behavior 

• Complex Preloads and Force Distributions 

• Fluid-Structure Interaction 

• Interaction with System Model 

Both the variations in the behavior and the complexity of the tire’s environment make 
it difficult for any one finite element analysis program to analyze all of the situations. 
The following discussion will outline each of the MSC product’s capabilities with 
emphasis on tire modeling and simulation. 



MSC/NASTRAN CAPABILITIES 


MSC/NASTRAN Version 68: 

• library of finite elements (scalar, 1-D, 2-D, 3-D, rigid, etc.) 

• h and p elements 

- can be used in same model 

• analysis capabilities: 

- linear static analysis 

- normal modes analysis 

- linear buckling analysis 

- transient and frequency response 

- nonlinear static and transient dynamics 

- complex eigenvalue analysis 

- steady state and transient heat transfer (linear and nonlinear) 

- superelement analysis 

- aeroelastic analysis 

- multidisciplinary design sensitivity and optimization 

• linear, nonlinear and hyperelastic materials 

• isotropic, orthotropic, and anisotropic materials 

• numerous performance enhancements 

• open architecture - interfaces to CAD and modeling systems 

• available on PCs, workstations, and supercomputers 

• fully documented and QA’d 

• on-line documentation 
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EXAMPLES OF TIRE SIMULATION USING MSC/NASTRAN 


For tire simulation MSC/NASTRAN ’s layered composite plate elements are widely 
used to model the ply lay-up. The slide line contact and gap elements simplify the contact 
and footprint characterizations (Figure 2). 




Figure 2. Sample Tire Model 



EXAMPLES OF TIRE SIMULATION USING MSC/NASTRAN 


Geometric and material nonlinearities are an important part of a tire simulation. The 
internal pressure preload in a tire affects its dynamic characteristics. The stiffening 
effects of a spinning tire can also be of importance. MSC/NASTRAN ’s geometric and 
material nonlinear capabilities can be used to include initial stress or preload conditions 
into a modal analysis. 

Once the initial internal pressure, wheel rim mounting loads and vertical load have 
been applied, the effects of applied lateral forces and cornering loads can be 
characterized. Figure 3 shows a cross section of a tire-wheel assembly subjected to a 
lateral force. 




Figure 3. Tire - Wheel Assembly Subjected to Lateral Force 
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EXAMPLES OF TIRE SIMULATION USING MSC/NASTRAN 


Figure 4 shows a close-up of the displaced shape of a tire cross section, first subjected 
to vertical load, and then subjected to a lateral force. The effects of contact area and 
friction force can be determined using these types of models. 



Figure 4. Displaced Shapes of Tire Cross Sections 



ADDITIONAL FEATURES OF MSC/NASTRAN 


MSC/NASTRAN has the following additional features for tire simulation: 




automatic time stepping in linear and nonlinear transient analyses, 
automatic load stepping in nonlinear statics, 
adaptive gap element 
slide line contact 

multidisciplinary design sensitivity and optimization 

- Simultaneous multiple loadings 

* static loads 

* modal results 

* dynamic results 

- design constraints 

* stress, displacement, natural frequency, etc. 

- physical properties as Design variables 

* plate thickness, composites ply thicknesses and orientation 
angles, etc. 

DMAP (Direct Matrix Abstraction Program) programming language 
superelement analysis 
cyclic-symmetry analysis 


138 



TIRE SIMULATION WITH MSC/ABAQUS 


MSC/ABAQUS 

• powerful nonlinear static and dynamic analysis 

• offered by MSC (agreement with Hibbitt, Karlsson, and Sorensen, Inc., HKS) 

• specific capabilities to address the nonlinear characteristics of tire problems. 

• nonlinear material models 

• slide line, gap and special contact elements 

• user routines for representing special materials 

• special purpose elements for unique applications. 

MSC/ABAQUS plate and solid Rebar elements (Figure 5) 

• initially developed for modeling reinforced concrete 

• easily adapted to tire problems 

• ply cords can have multiple levels (plies), each with different properties and 
orientations. 

• material properties can be nonlinear 

• can be combined with other nonlinear (i. e., hyperelastic) materials that form 
the tire tread stock. 




Figure 5. MSC/ABAQUS Rebar Elements 


A number of MSC/ABAQUS’s nonlinear analysis capabilities are included in 
MSC/PATRAN Advanced Finite Element Analysis (AFE A). AFEA is a powerful 
tool combining a subset of MSC/ABAQUS capabilities with MSC/PATRAN as the 
integrated graphical pre- and postprocessor. 
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MSC/DYTRAN FOR TIRE SIMULATION 


MSC/DYTRAN 

• highly nonlinear short duration transient dynamic analysis 

• library of finite elements 

- scalar 

- 1-D 

- 2-D 
-3-D 
-rigid 

• Lagrange processor for structural applications 

• Euler processor for modeling fluids 

• special elements for fluid/structure coupling. 

• explicit integration 

- short duration events (such as impact) 

• complex general contact 

• special elements 

• nonlinear material models 

• Input similar to MSC/NASTRAN 

- can use MSC/PATRAN and SDRC I-DEAS 

• User routines for special conditions 

• Link to MSC/NASTRAN for static or quasi-static initial conditions 
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MSC/DYTRAN FOR TIRE HYDROPLANING SIMULATION 


Tire hydroplaning simulation (Figure 6): 

• rotating tire 

• specified velocity 

• lowered into a fluid 

• displacement condition - added to simulate a turning maneuver. 

• different road surfaces were considered. 

• coarse mesh used for demonstration purposes 

• valuable insight into behavior 

• production level analysis requires a finer mesh. 



Figure 6. Tire Model 


Tire model - Lagrangian elements 

• Shell elements for tire sidewalls 

• Solid elements for tread 
Wheel - assumed rigid 

• Set of rigid elements 
Water model - Eulerian elements 
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MSC/DYTRAN FOR TIRE HYDROPLANING SIMULATION 


The fluid region of the tire model is a combination of elements containing water and 
void elements. A cut-away view of the model is shown in Figure 7. As the tire moves 
through the water, the water will rise around the tire into the void elements of the 
Eulerian mesh. 



Figure 7. Cutaway View of Tire Fluid - Structure Model 
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MSC/DYTRAN FOR TIRE HYDROPLANING SIMULATION 


As the simulation progresses, the water rises and lifts the tire off the road surface. 
MSC/DYTRAN ’s Euler processor enables the engineer to track the motion and pressure 
distributions of the fluid around the tire. In the simulation, water moves through the 
Euler mesh into areas of the mesh that were initially a void. Figure 8 shows sample fluid 
pressure contours superimposed on the fluid displaced shape. 



Side View End View 


Figure 8. Sample Fluid Displacement and Pressure Contour Plots 
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MSC/DYTRAN FOR TIRE HYDROPLANING SIMULATION 


The effect of the fluid forces lifting the tire off the road surface can be studied by 
plotting the time variation in the vertical deflection of a node at the center of the wheel 
(hub) as shown in Figure 9. 


1-02 



Figure 9. Time History Plot of Wheel Hub Vertical Deflection 
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SUMMARY 


Tires are complex structures that experience a wide range of loadings under varied 
conditions during their service life. The MacNeal-Schwendler Corporation has many finite 
element based structural analysis products available (MSC/NASTRAN, MSC/ABAQUS, 
MSC/PATRAN AFEA, and MSC/DYTRAN) to address the many complex aspects of tire 
behavior. Many of the capabilities of these products can be uniquely applied to tire 
simulation. The availability of User routine customization is also possible to address 
special materials, special loads and special boundary conditions. 
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and 
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201 Johnson Road 
Houston, PA 15342-1300 


INTRODUCTION 

This presentation will provide a summary of the capabilities in the ANSYS program that 
relate to the computational modeling of tires. The power and the difficulties associated with 
modeling nearly incompressible rubber-like materials using hyperelastic constitutive relationships 
will be highlighted from a developer’s point of view. The topics covered in this presentation will 
include: a hyperelastic material constitutive model for rubber-like materials, a general overview of 
contact-friction capabilities, and the acoustic fluid-structure interaction problem for noise 
prediction. Brief theoretical development and example problems will be presented for each topic. 


149 



POSSIBLE FACTORS AFFECTING TIRE MODELING 


Some of the factors that directly and indirectly affect the design and analysis of tires are listed 
below. The ANSYS program has extensive capabilities in each of these areas. However, in this 
presentation, attention will be focussed on the material modeling of rubber and rubber-like 
materials. Specifically a hyperelastic material model will be studied. The ability of ANSYS to solve 
contact problems with friction will also be briefly summarized. Finally, since noise control is of 
particular interest to tire manufacturers, a brief overview of the ANSYS capabilities of solving 
acoustic fluid-structure interaction problems is given. 


• Modeling of rubber or rubber-like materials - hyperelasticity <= 

• Contact-Friction Capability 

• Composite Materials 

• Coupled Acoustic Analysis 

• Coupled Thermal Analysis 

• Other Design Tools, e.g.. Design Optimization, APDL, etc. 

Figure 2 
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MOONEY-RIVLIN HYPERELASTIC MATERIAL MODEL 


The ANSYS program offers the Mooney-Rivlin hyperelastic material model to characterize 
rubber and rubber-like materials. This material model is suitable for capturing incompressible and 
nearly-incompressible behavior of rubber-like materials. Up to nine material constants can be used 
to define the strain energy density stored in a material. These constants are derived from laboratory 
tests on rubber specimens. A procedure for automatically determining these constants from 
experimental data will be discussed later. The Mooney-Rivlin hyperelastic material law is available 
with the ANSYS elements HYPER56 (4-noded quadrilateral), HYPER58 (8-noded brick) and 
HYPER74 (8-noded quadrilateral). 


• Mooney-Rivlin is a material option applicable to incompressible rubber or rubber-like 
materials. Two, five, or nine material constants are used to describe the material behavior. 


If the Mooney-Rivlin constants for a material to be analyzed are known, the hyperplastic 
constants are input into ANSYS with the TB and TBDATA commands. If the 
Mooney-Rivlin constants are unknown for a rubber material, they can be derived from 
laboratory test data and directly sent to ANSYS using the *MOONEY command. The 
procedure for determining the constants from test data will be discussed later. 


• The Mooney-Rivlin option is available with the incompressible elements HYPER56, 58, 
and 74. Up to six temperature-dependent sets of Mooney-Rivlin constants may be 
implemented. 


Figure 3 
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MOONEY-RIVLIN STRAIN ENERGY DENSITY FUNCTIONS 


The strain energy density function of Mooney-Rivlin materials is normally expressed as a 
function of the invariants, I; , of the right Cauchy-Green tensor, C, and a number of constants. In 
order to separate the distortional and dilatational parts of deformation, a set of reduced invariants 
(defined in the next transparency), Jj, are used in place of the original invariants, Ij. The constant x is 
the bulk modulus of the material and is specified through the Poisson’s ratio v. The value of x 
controls the degree of incompressibility of the materials. At v = 0.5, x is infinity and the material is 
fully incompressible. A value of v, very close to 0.5, viz., 0.4999, can be used in the program to 
model nearly incompressible rubber-like materials. Up to nine constants are allowed and up to six 
temperature-dependent sets of constants may be input. This representation of the strain energy 
density function can handle up to 400% strains as long as the element distortions are of acceptable 
limits. 


For the elements HYPER56, HYPER58, and HYPER74, the two, five and nine parameter 
Mooney-Rivlin models are available. The strain energy density function for these elements is given 
in polynomial form by: 

N 

W = X a kL (J 1 - 3)k (J 2 - 3 > L + 1 / 2lC ( J 3 - D 2 
K + L= 1 

where a^ = the constants of the nine-parameter cubic Mooney-Rivlin 

relationship 

x = ^ = bulk modulus 

(l-2v) 

Setting N=3 in the polynomial Mooney-Rivlin strain energy density expression, we obtain 
the nine parameter Mooney-Rivlin model: 

W = a 10 (Ji - 3) + aoi (J 2 - 3) + a 20 (h - 3) 2 + an (Jj - 3)(J 2 -3) 

+ ao 2 (J 2 - 3) 2 + a3o (Ji - 3) 3 + a 2 i (Ji - 3) 2 (J 2 - 3) 

+ a 12 (J 1 -3)(J 2 -3)2 + a 03 (J 2 -3)3 + l/ 2 x(J 3 -1)2 

Note that analogous two and five parameter models are obtained by respectively setting N= 1 
and N=2 in the polynomial expression. The last term in the model represents the hydrostatic 
(volumetric) work. 


Figure 4 
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REDUCED OR VOLUME-FILTERED INVARIANTS 


Figure 5 defines the reduced invariants, Jj, in terms of the original invariants, I;, of the right 
Cauchy-Green tensor. Note that the invariants are being scaled by the determinant of the 
deformation gradient, det (F) = 13, which also represents the change of volume of the material. In the 
hyperelastic formulation, ANSYS uses the total Lagrangian approach in which the original 
coordinate system is always used as the frame of reference. The energy conjugate in this case is: 
(C:S), where S is the second Piola-Kirchhoff stress. The strain C and the stress S are transformed 
internally and the results are presented in terms of Hencky strain and true stress. 


j, = i,i 3 ->/ 3 

J 2 = I 2 I 3 ~ 2/3 
J 3 = Det (F) = I 3 ‘/ 2 


where: 


F 

If *3 


Deformation Tensor 

Invariants of right Cauchy-Green tensor C = F T F 


Figure 5 
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INCOMPRESSIBILITY 


Rubber-like materials undergo large deformations with small changes in volume. The 
incompressible or nearly incompressible behavior of hyperelastic materials becomes more and more 
pronounced as Poisson’s ratio, v, approaches 0.5. The incompressibility condition must be 
accounted for in the element formulation so that the element exhibits proper behavior under loading. 
Correct treatment of the incompressibility condition in the element formulation not only prevents 
the element from locking under certain loading and boundary conditions, but also enhances the 
convergence characteristics of the element. 


• Rubber-like materials undergo volume-preserving deformations. 


• If Poisson’s ratio v is close to 0.5 for hyperelastic materials, then this incompressibility 
condition must be accounted for in the element formulation. 


Figure 6 



INCOMPRESSIBILITY EQUATION 


The hyperelastic elements in the ANS YS program use extra pressure degrees of freedom in 
an attempt to satisfy the incompressibility condition. An augmented strain energy density function 
with the pressure appearing explicitly is used instead of the original energy function. The constraint 
condition, which is shown below, is satisfied over an element in an average sense. Since J3 

A 

represents the volume change of an element derived from its displacements, and J 3 is the volume 
change derived from its computed pressure, it follows that the integrated volume change over an 
element would be maintained approximately at zero. 


The hyperelastic elements use pressure DOFs by augmenting the potential energy equation 

with: 

| {N p } T (J 3 - J 3 ) dV = 0 

{N p } are the pressure shape functions 

J3 is the volume change from the displacements 

A 

J 3 is the volume change from the computed pressure 
J 3 = det (F) 



x = bulk modulus and p = computed pressure 


Figure 7 



ELEMENT STIFFNESS MATRIX 


The element stiffness matrix equation for an element is shown below. The pressure degrees 
of freedom (dof), which are not assigned to any node, are eliminated at the element level. Therefore, 
the second equation, which is in terms of the displacements only, is sent out for assembly and 
solution. The meanings of the symbols in Figure 8 are: K uu = stiffness sub-matrix corresponding to 
the displacement dof, K pp = stiffness sub-matrix corresponding to the pressure dof, K pu = stiffness 
sub-matrix corresponding to the coupled displacement and pressure dof, Au = incremental 
displacement, Ap = incremental pressure, F a PP = applied external loads, F nr u = Newton-Raphson 
restoring force corresponding to the displacement dof, F nr p = Newton-Raphson restoring force 

A 

corresponding to the pressure dof, K = final element stiffness matrix modified by the condensation 

a nr 

of Ap, and F = final Newton-Raphson restoring force corresponding to the displacement dof 
modified by the condensation of Ap. 


This leads to an equation of the form 



In all cases, the pressures are interpolated as “nodeless” variables, allowing the equation to 
be condensed to displacement DOF only: 

a nr 

[K] {Au} = { F app } - {F } 


Figure 8 


156 



A FAMILY OF ANSYS HYPERELASTIC ELEMENTS 


As mentioned earlier, ANSYS offers a family of elements that use the Mooney-Rivlin 
hyperelastic constitutive model. HYPER56 is a 4— noded quadrilateral element with one pressure 
dof, HYPER74 is an 8-noded quadrilateral element with three pressure dof, and HYPER58 is an 
8-noded brick element with one pressure dof. The users are not concerned with the “nodeless” 
pressure dof which are eliminated at the element level. Furthermore, the 8-noded quadrilateral is 
actually a 9-noded Lagrangian element. The mid-side node is eliminated at the element level. Note 
that the pressure is interpolated using shape functions that are one order lower than the displacement 
shape functions. 



u + p 
4+ 1 


8 + 3 



• - displacement nodes 


o - pressure nodes 


Figure 9 






MOONEY-RIVLIN CONSTANT DETERMINATION 


The hyperelastic constitutive models are empirical in nature. The Mooney-Rivlin model is 
given in polynomial form in terms of the strain invariants and a number of unknown material 
constants. One of the difficulties of using the Mooney-Rivlin model is the uncertainty associated 
with ascertaining the values of the unknown constants. These material constants cannot be chosen 
arbitrarily; they should satisfy the constraints listed in Figure 10. If one or more of the constraints are 
violated, the program issues warning messages. The constraints are usually violated when 
insufficient experimental data was used to derive the constants. As long as the analysis is within the 
strain range of the test data and undergoes the same mode(s) of deformation (e.g., uniaxial tension), 
then the constraints may be violated without affecting the solution accuracy. If the constraints are 
violated and the model experiences modes of deformation which there is no experimental data for, 
then convergence difficulties may be encountered. In this case, the results should be carefully 
validated. 


The material constants are not arbitrary, but must satisfy certain constraints in order to 
produce valid behavior in all modes of deformation. These constraints are imposed in order to 
maintain: 

• Positive definiteness of the strain energy density function 

• Stress is always positive (or negative) at the extremes of uniaxial deformation 

• Stress is a continuous function of deformation 

If any constraint is not satisfied, the program will issue a warning message. 


Figure 10 
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PROCEDURE FOR DETERMINING MOONEY-RTVLIN CONSTANTS 


Laboratory tests are performed on the representative rubber test specimens to obtain 
experimental stress-strain data. Uniaxial tension and compression tests are the most common, but 
equibiaxial tension, equibiaxial compression, and shear tests are also regularly performed. The 
ANSYS program takes input engineering stress-strain data and performs a least-square 
curve-fitting procedure to determine the Mooney-Rivlin constants of the hyperelastic strain energy 
density function. Experimental data from more than one test may be used to determine the 
Mooney-Rivlin constants. The user has the option of determining the constants without performing 
analysis or combining the constant determination with detailed analysis. The program automatically 
stores the Mooney-Rivlin constants in the ANSYS database, and a user specified data array, and 
writes them to an output file in typical ANSYS data table format. 




The Mooney-Rivlin constants for hyperelastic materials are not generally available in the 
open literature. Therefore, most sets of Mooney-Rivlin constants need to be derived from 
experimental stress-strain data. 

For hyperelastic materials, simple deformation test data can be used to accurately 
characterize the Mooney-Rivlin constants. The constants are determined by ANSYS using 
the *MOONEYcommand. 

Test data may come from one or more of the following six tests: 

— Uniaxial Tension - Equibiaxial Tension 

- Uniaxial Compression - Equibiaxial Compression 

- Planar Tension (Shear) - Planar Compression (Shear) 


Figure 11 
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EXPERIMENTAL STRESS-STRAIN DATA 


Figure 12 shows two sets of experimental data, one from a uniaxial tension test and the other 
from a uniaxial compression test. The data was obtained from experimental tests conducted on 
rubber specimens taken from the same stock. The stress-strain values are expressed in engineering 
form for practicality. The ANSYS program can determine Mooney-Rivlin hyperelastic material 
constants from experimental data obtained from any combination of the three independent modes of 
deformation: uniaxial tension, equibiaxial tension and planar tension (shear). By superimposing 
hydrostatic pressure to the loaded specimen, uniaxial compression tests can be shown to be 
equivalent to equibiaxial tension tests. Similarly, uniaxial tension data can be equated to equibiaxial 
compression data and pure shear (planar tension) data can be equated to planar compression data. It 
is noteworthy to mention that Mooney-Rivlin constants cannot be independently determined from 
shear test data alone; equations from another test, viz., uniaxial tension or equibiaxial tension, must 
also be considered with the shear test equations. 


COMPRESSION DATA 

TENSION DATA 

strain 

stress* 

strain 

stress* 

-.4968 

-1345.74 

.0024 

5.13238 

-.4708 

-1174.01 

.0142 

15.39714 

-.4430 

-1020.92 

.0317 

28.22864 

-.4162 

-897.182 

.0507 

42.34296 


... 

.0685 

56.45728 

-.0415 

-50.0036 

.0840 

71.85739 

-.0356 

-41.2488 

.1012 

83.40508 

-.0297 

-32.7781 

.1184 

94.95266 

-.0247 

-25.4332 

.1353 

105.21460 

-.0197 

-18.6461 

.1528 

115.48130 

-.0155 

-13.5637 



-.0121 

-9.6072 



-.0088 

-6.2189 



-.0054 

-3.3988 



-.0037 

-1.7047 



-.0020 

-.2841 



* stress in psi 


Figure 12 
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PLOTS OF EXPERIMENTAL AND FITTED CURVE 


After the constants have been determined, the user should plot the experimental data curve 
and the computed curve on the same frame in order to view the quality of the fit. Figure 13 shows 
such a plot. Since the two curves are coincident, it appears that the fit in this case was adequate. As a 
rule of thumb, the number of data points should be at least twice the number of constants requested in 
order to obtain a reasonable fit. For the present example, there were 42 data points for the 
compression test and 10 data points for the tension test. A nine-parameter curve was fitted. The 
coefficient of determination, which reflects how well the curve fits the input data, is also printed in 
the output file. The closer the value to unity, the better the fit. However, a value closer to unity does 
not necessarily imply that the curve is “good”, only that it goes through the input points fairly well. 
The coefficient of determination for the present example was 0.9998. Another measure of the 
quality of the fit, the root-mean-square error, is also calculated and output. The root-mean-square 
error, expressed as a percentage, should be “close” to zero. Its value for this example was 8.53%. 



Figure 13 



EVALUATING MOONEY-RIVLIN CONSTANTS 


In graphing the calculated stress-strain curve in the ANSYS program, it is possible to extend 
the displayed curve into regions that were not defined by the experimental data. Graphing the curves 
over such an extended range can help the user qualitatively understand the model’s behavior if its 
response ever happens to exceed the range of the experimental strain. It should be noted that a 
uniaxial tension curve should be graphed only in regions of positive strain, and a uniaxial 
compression curve only in regions of compressive strain. A good practice usually requires that the 
test data represent all modes of deformation and ranges of response experienced by the model. 


• In order to evaluate the calculated constants between the data points or to analyze the 
stress-strain behavior outside the range of test data, the *Evaluate command should be used. 

• To establish if there are any spurious inflections in the calculated stress-strain curve between 
data points, the determined Mooney-Rivlin constants should be evaluated in at least 1 0 times 
the number of data points. 

• When extrapolation from the data points in a hyperelastic model is necessary, the *Evaluate 
command should always be used to verify expected material behavior. 


Figure 14 
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PLOT OF EVALUATED STRESS-STRAIN CURVE 


Figure 15 shows a portion of the calculated stress-strain curve that goes beyond the range of 
input data. The range of input data for this example lies in the middle part of the graph which is 
extended in the tension region. In the tension region, beyond the range of input data, the material is 
expected to exhibit a brief plateau followed by a stiff strain-hardening behavior. This knowledge 
may become important in the interpretation of results of an analysis in which parts of the model 
experienced tensile strains in this range. 



Figure 15 
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INCOMPRESSIBILITY: LOCAL AND GLOBAL B-B CONDITIONS 


For a mixed formulation with displacement and pressure dofs, finite elements can be 
developed that satisfy the so-called Brezzi-Babuska (B-B) inf-sup condition. The finite elements 
that satisfy the Brezzi-Babuska condition are not expected to lock at the element level. However, the 
satisfaction of the B-B condition at the element level does not always guarantee a lock-free behavior 
at the global level. A problem with a certain set of boundary and loading conditions can always be 
created under which the problem will exhibit a locking behavior. Figure 16 gives guidelines for a 
proper mix of displacement and pressure dof in a model. Also, the analyst should watch out for the 
occurrence of actual buckling in his problem. 


Guidelines: 

two-dimensional Problems - No. of free uj > 2pj 
three-dimensional Problems - No. of free u; > 3pj 

Also, 

•sp Watch for real buckling 

Figure 16 
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USER-DEFINED MATERIAL LAWS 


Tire manufacturers may have their own material laws for describing rubber-type material 
behavior which are proprietary in nature. ANSYS provides user-programmable features, e.g., a 
hyperelastic user-routine, which allows direct input of personalized material constitutive laws into 
the program. The advantage of this procedure is that the extensive pre-processing, solution, and 
post-processing capabilities of ANSYS can be utilized without revealing a company’s proprietary 
material constitutive model. Because of competition in the tire industry, the user-programmable 
features of ANSYS could prove very useful to the concerned companies. 


If there is a need to simulate a material behavior that is not a program option, then the user 
may want to consider programming the material law himself. 

A number of user-programmable features (UPFs) are offered for characterizing different 
types of material behavior. 

The available user material routines include: 

• Plasticity 

• Hyperelasticity <= 

• Creep 

• Swelling 


Figure 17 
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USER-PROGRAMMABLE ROUTINE FOR HYPERELASTIC MATERIAL 


Figure 18 shows the structure of the hyperelastic user-routine. The routine provides the 
invariants of the right Cauchy-Green tensor, Poisson’s ratio and problem dimensionality as the 
input. The user is required to compute the strain energy density, the pressure, the first and second 
derivatives of the strain energy density function and the pressure with respect to the invariants, and 
the bulk modulus. This structure not only provides the mechanism for writing an invariant-based 
hyperelasticity material law, but also a principal-stretch-based hyperelasticity material law. To 
write code for a constitutive law of the latter type: (1) compute the principal stretches using the 
invariants, (2) compute the first and second derivatives of the strain energy density function (which 
is now a function of the principal stretches) with respect to the principal stretches, and (3) transform 
the derivatives with respect to the principal stretches into those with respect to the invariants using a 
chain rule of differentiation. 


Input quantities provided by the program: 


1. Pl,l2,I 3 ] 


2. Poisson’s Ratio 


3. Problem Dimensionality 
[ 2 = two-dimensional, 3 = three-dimensional ] 


User responsibility is to compute: 


1. W, 


aw a^w 
di ’ ai 2 


P ap a 2 P 
’ si ’ ai 2 


3. K = Bulk Modulus 


Figure 18 
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SNAP-THROUGH OF A HYPERELASTIC BUTTON 


As previously described, the ANSYS program allows hyperelastic materials to undergo 
realistic, large recoverable deformations. The geometry in Figure 19 depicts the buckling of a soft 
silicon rubber button that is used in a mouse to activate electrical impulses. The button, which was 
modelled using axisymmetric HYPER74 elements, actually ’snaps-through’ as it moves from its 
initial to final configurations. The deformed shape demonstrates the hyperelastic element’s ability 
to undergo extremely large strains, as is shown in the deepening folds which physically occur as the 
button displaces. 



(a) Initial Geometry 



(b) Deformed shape with stress contour 


Figure 19 
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PRESSURIZED AXISYMMETRIC O-RING 


In order to demonstrate how hyperelastic materials can be used in conjunction with contact 
elements in ANSYS, a hyperelastic O-ring seated in a square groove is analyzed. The O-ring of 
solid circular cross-section is first forced to conform to the shape of a channel by installing an inner 
retaining ring against the O-ring. Pressure is then applied to the top surface of the O-ring. In the 
model, the O-ring is meshed with axisymmetric HYPER56 elements and the rectangular channel 
and moving retaining ring consist of two-dimensional spar (LINK1) elements. Two-dimensional 
point-to-surface contact elements (CONTAC48) are used to define contact between the O-ring and 
the spar elements. As expected, the analysis of the O-ring finds that the maximum stress occurs at 
the corners of the channels. 



(a) Initial geometry (b) Deformed shape with hydrostatic stresses 



Figure 20 
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GENERAL CONTACT CAPABILITIES 


In the ANSYS program, node-to-node (nodal) contact is represented by elements 
CONTAC12, 52, and COMBINE40 while surface-to-surface (general) contact is represented by 
elements CONTAC26, 48, and 49. For nodal contact, the final points of contact must be known 
beforehand. In general surface contact, contact areas do not have to be known a priori and large 
sliding motions and deformations are permitted. In fact, almost all nonlinear options - plasticity, 
large strain, hyperelasticity, etc. - can be used in conjunction with general surface contact. Since the 
surface contact elements are defined using any pre-existing solid element (including shells) mesh, 
simple and complex planar, axisymmetric, and three dimensional configurations can be modelled 
with contact. In ANSYS, CONTAC26 and 48 are used to model two-dimensional and 
axisymmetric surface contact while CONTAC49 is used to define contact between three 
dimensional elements. Surface contact with midside-node elements is also allowable in ANSYS 
when additional steps are taken to alleviate non-uniform stiffness along mid-sided surface nodes. 
All three general surface contact elements can be used in static, dynamic, and thermal-structural 
analyses. 


• In ANSYS, both point-to-point and general surface-to-surface contact elements are 
available. In all cases, friction may be included at the contact interface. 

• The points or surfaces may be either rigid or deformable so that the component may come 
into contact with a rigid element or itself. Note that point-to-surface contact is a special case 
of surface-to-surface contact. 

• Large deformations, plasticity, and other nonlinear effects can be used in conjunction with 
frictional contact. 

• Planar, axisymmetric, and three-dimensional configurations (including shells) are allowed 
in contact analyses. 

• Contact surface elements are defined using pre-existing element meshes so that complex 
surfaces can be easily represented. 

• Static, dynamic, and coupled thermal-structural contact analysis options are available. 


Figure 21 



CONTACT SURFACE KINEMATICS 


When using general surface contact elements, ANSYS performs a multi-level search for 
contact to increase program efficiency. Utilizing pinball regions whose radii are defined by the 
length of each contact plane, four levels of contact status are established. When a contact node is far 
from the contact plane well outside the pinball region, STAT=4. For such far field status, the contact 
algorithm is simple and utilizes minimal CPU time. For a node near the contact plane inside the 
pinball region, but out of contact, STAT=3. In this case, the contact algorithm is more complex and 
uses a modest amount of CPU. Once contact has been established between the nodes of two surfaces, 
a STAT value of 1 or 2 is given to the contact element. STAT=1 corresponds to a node being in 
contact and sliding while STAT=2 signifies a node being in contact and sticking. When contact 
occurs, the contact algorithm is very complex and uses a large amount of CP time. 


The kinematics of contact deal with the position and motion of the contact nodes relative to 

the contact planes. In surface-to-surface contact analyses: 

• A multi-level search procedure is performed in ANSYS to establish if contact has occured. 
In this procedure, a distinction between near-field and far-field states of contact is made 
using the pinball concept. 

• Contact voids and overlaps are eliminated since ANSYS generates pseudo-contact elements 
around each contact plane in such a manner that each contact node is uniquely defined with 
respect to a target surface. 

• Both symmetric and asymmetric contact analyses options are available. 


Figure 22 



SUMMARY OF CONTACT SURFACE ANALYSIS WITH FRICTION 


The tangential sliding behavior of two contacting bodies in relative motion may be 
frictionless or may involve friction. Frictionless behavior allows the two bodies to slide relative to 
one another without any resistance. Including friction, however, causes shear forces to develop 
between the contact regions of the two bodies. ANSYS is capable of defining contact with and 
without friction. Frictionless contact is suitable in models whose components have low friction 
coefficients and develop minimal tangential forces. Two dry friction models are available in 
ANSYS for surface-to-surface contact: rigid and elastic Coulomb friction. Rigid Coulomb friction 
allows only a sliding friction condition where the contact areas do not consider sticking. Elastic 
Coulomb friction allows both sticking and sliding conditions. The sticking zone is treated as an 
elastic region in which the components in contact will return to their original geometry when the 
tangential force is removed. At forces above static friction, the components will slide inelastically 
relative to each other. Elastic Coulomb friction is valid for most engineering materials in contact. A 
user defined friction option is also available for modeling friction behavior which is dependent on 
velocity. This option is most commonly used for cases where fluid-film effects are important. 
Fluid-film effects are characterized by a change in the tangential friction force with an increase in 
contact surface velocity. 


In ANSYS, four different friction models are available for general surface-to-surface 
contact: 

i) frictionless contact behavior 

ii) dry elastic Coulomb friction 

iii) dry rigid Coulomb friction 

iv) a user-friction capability that allows the influence of tangential velocity on the friction law. 

Figure 23 
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ANALYSIS OF A HEADING PROCESS 


To illustrate how contact elements can be used with other nonlinear analysis techniques, the 
heading process of a rivet was analyzed using ANSYS. In the analysis, an aluminum workpiece of 
circular cross section was press-fit into a steel grip die using a hemispherical heading tool. The 
workpiece, grip die, and heading tool were all modelled using axisymmetric PLANE42 elements. 
CONTAC48 elements, with the rigid Coulomb friction option, were used to simulate the contact 
between the workpiece, heading tool, and grip die. Since the workpiece undergoes large 
nonrecoverable deformations, plasticity was included in the model. As demonstrated by the 
deformed shape, the results of the analysis accurately predict that the maximum stress occurs along 
the comer of the grip die which the workpiece rolls over. 



Figure 24 
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TAPER SHAFT HEAT UP CONTENT 


To demonstrate the coupled thermal-structural capabilities of ANSYS, the expansion of a 
heated tapered shaft is analyzed. A coupled thermal-structural analysis uses coupled field solids 
PLANE 13, SOLID5, or SOLID98 elements. These elements have both thermal and structural 
degrees of freedom. Heat flow occurs in these elements when contact is established and the 
CONTAC48 or 49 elements have thermal-structural degrees of freedom activated. In the tapered 
shaft analysis which uses SOLID98 and CONTAC49 elements, both temperature and stress contours 
are plotted to demonstrate the thermal-structural coupling. 


FLUID-STRUCTURE COUPLING 


Noise prediction and control are also important considerations for the tire manufacturers. 
The ANS YS program has extensive capabilities for predicting pressure levels generated by vibrating 
structures or other sound sources. In order to solve for the pressure distribution generated by 
vibrating structures, the structural dynamic equations must be coupled to the acoustical fluid flow 
equations through the use of the momentum balance equation at the fluid-structure interface. 
Depending upon the nature of the excitation source, the problem can be solved in the time domain or 
in the frequency domain. The ANSYS program has separate acoustic absorption elements to model 
acoustic problems whose domain extends to infinity. If acoustic absorption elements are used, the 
domain is truncated at a finite distance and the absorption elements are placed at the truncated 
boundary. 


• Dynamic Acoustic Fluid Equation 

[K f ] {p} + [M f ] {p} = 


4ds 

dn 


Structural Dynamic Equation 


[M s ] {ii} + [C s ] {u} + [K s ] {u} = 
• Fluid-Structure Interface Condition 


Nj p ds + {F s } 


dp _ „ d 2 u 

dn e ° at 2 


• Combining the Fluid-Structure Equations 
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where: 


R 

= coupling matrix 



Cf 

= Fluid damping matrix at the boundary surface to account for sound 
absorption 



00 

= fluid density 
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PASSENGER CAR MODEL 


A finite element model of an automobile body with seats is shown. The interior acoustic 
enclosure is modelled using three-dimensional acoustic fluid elements (FLUID30) and the car body 
is modelled using shell elements (SHELL63). The acoustic elements in touch with the enclosure 
walls are coupled with the structural shell element. The velocity of sound travelling through the air 
in the automobile enclosure is taken to be 13200 inches per second, and the density of air 0.0001958 
pounds-second**2 per inch. Two vibrating panels, one on the driver’s side and the other on the 
passenger’s side, simulate loudspeakers of the automobile audio-system. A harmonic analysis is 
performed in which the loudspeakers are excited at a frequency of 40 hertz. The objective is to 
monitor the pressure levels at the locations of the riders’s ears. A contour plot of the absolute 
pressure on the passenger side of the automobile compartment is shown. The analysis can be carried 
out at different operating audio frequency levels. The present analysis did not consider damping on 
the enclosure walls. Empirical damping values could be easily added to make the analysis more 
realistic. 




Location of monitoring 
the pressure response 


Back 


Seat 


Location of the 


Front 

Seat 


vibrating panel 


(a) Perspective View 


(b) Pressure Plot on 
Passenger Side 
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